Bamboo mesocosms experiment
11 April 2025
1 Setup
1.1 Libs
#install.packages("remotes")
#install.packages('TMB',source=TRUE) #didn't work on home laptop
#remotes::install_github("glmmTMB/glmmTMB/glmmTMB")
library("glmmTMB")## Warning in check_dep_version(dep_pkg = "TMB"): package version mismatch:
## glmmTMB was built with TMB package version 1.9.15
## Current TMB package version is 1.9.16
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)#install.packages("BiocManager")
#BiocManager::install("phyloseq")
library("phyloseq")
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union##
## Attaching package: 'reshape'## The following object is masked from 'package:dplyr':
##
## rename## Loading required package: stats4##
## Attaching package: 'bbmle'## The following object is masked from 'package:dplyr':
##
## slice## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:reshape':
##
## rename, round_any## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize1.2 Re-read data
setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Bamboo_mesos/Bamboo_mesos/04.mixed_models/")
#ps.clean <- readRDS("../01.process_asvs/reusum23.ps.lulu.clean.rds")
ps.clean <- readRDS("../01.process_asvs/bamboomesos_ps.lulu.clean.decontam.rds")
#ps.clean <- readRDS("../01.process_asvs/ps.all.clean.100.rds")
ps.clean #514 taxa, 169 samples## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 514 taxa and 169 samples ]
## sample_data() Sample Data: [ 169 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 514 taxa by 10 taxonomic ranks ]ps.exp <- subset_samples(ps.clean,Mesocosm_type=="Experiment"&Microbe_treatment!="AHK"&Exp_day!=20&Microbe_treatment!="ECO"&Microbe_treatment!="ALL")
ps.exp #60 samples## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 514 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 514 taxa by 10 taxonomic ranks ]##including ECO & ALL to compare
ps.exp.all <- subset_samples(ps.clean,Mesocosm_type=="Experiment"&Microbe_treatment!="AHK"&Exp_day!=20)
ps.exp.all #100 samples## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 514 taxa and 100 samples ]
## sample_data() Sample Data: [ 100 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 514 taxa by 10 taxonomic ranks ]## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 514 taxa and 30 samples ]
## sample_data() Sample Data: [ 30 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 514 taxa by 10 taxonomic ranks ]##plus eco & all
ps.exp.all.wtr <- subset_samples(ps.exp.all,Sample_type=="Water")
ps.exp.all.wtr #50 samples## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 514 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 514 taxa by 10 taxonomic ranks ]## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 514 taxa and 30 samples ]
## sample_data() Sample Data: [ 30 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 514 taxa by 10 taxonomic ranks ]##plus eco & all
ps.exp.all.lar <- subset_samples(ps.exp.all,Sample_type=="Larvae")
ps.exp.all.lar #50 samples## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 514 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 514 taxa by 10 taxonomic ranks ]2 Data formatting
2.1 Trim, thoroughly
Ran once then saved
##6 samples out of 50 = 12%
##6 samples out of 30 = 20%
##water
ps.exp.wtr.trim <- filter_taxa(ps.exp.wtr, function(x) sum(x > 1) > (0.20*length(x)), TRUE)
ps.exp.wtr.trim ##62 taxa## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 62 taxa and 30 samples ]
## sample_data() Sample Data: [ 30 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 62 taxa by 10 taxonomic ranks ]## WB9 WL13 WB7 WL6 WL23 WL9
## 77861 72372 68692 63574 53741 51123##with all and eco
ps.exp.all.wtr.trim <- filter_taxa(ps.exp.all.wtr, function(x) sum(x > 1) > (0.12*length(x)), TRUE)
ps.exp.all.wtr.trim #69 taxa## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 69 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 69 taxa by 10 taxonomic ranks ]##larvae
ps.exp.lar.trim <- filter_taxa(ps.exp.lar, function(x) sum(x > 1) > (0.20*length(x)), TRUE)
ps.exp.lar.trim #44 taxa ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 44 taxa and 30 samples ]
## sample_data() Sample Data: [ 30 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 44 taxa by 10 taxonomic ranks ]##with all and eco
ps.exp.all.lar.trim <- filter_taxa(ps.exp.all.lar, function(x) sum(x > 1) > (0.12*length(x)), TRUE)
ps.exp.all.lar.trim #54 taxa## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 54 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 54 taxa by 10 taxonomic ranks ]## LB8 LB9 LL12 LL7 LB11 LB12
## 36318 34226 33371 31971 18211 13437##proportion of zeroes
##water
otu.exp.wtr <- data.frame(ps.exp.wtr@otu_table)
sum(otu.exp.wtr == 0) / (ncol(otu.exp.wtr) * nrow(otu.exp.wtr)) * 100## [1] 93.19066#93.19066%
otu.exp.wtr.trim <- data.frame(ps.exp.wtr.trim@otu_table)
sum(otu.exp.wtr.trim == 0) / (ncol(otu.exp.wtr.trim) * nrow(otu.exp.wtr.trim)) * 100## [1] 54.94624#54.94624%
##water - with all
otu.exp.all.wtr <- data.frame(ps.exp.all.wtr@otu_table)
sum(otu.exp.all.wtr == 0) / (ncol(otu.exp.all.wtr) * nrow(otu.exp.all.wtr)) * 100## [1] 94.27237#94.27237%
otu.exp.all.wtr.trim <- data.frame(ps.exp.all.wtr.trim@otu_table)
sum(otu.exp.all.wtr.trim == 0) / (ncol(otu.exp.all.wtr.trim) * nrow(otu.exp.all.wtr.trim)) * 100## [1] 64.11594#64.11594%
##save if ready
#saveRDS(ps.exp.wtr.trim,file="ps.exp.wtr.trim.rds")
#saveRDS(ps.exp.all.wtr.trim,file="ps.exp.all.wtr.trim.rds")
##larvae
otu.exp.lar <- data.frame(ps.exp.lar@otu_table)
sum(otu.exp.lar == 0) / (ncol(otu.exp.lar) * nrow(otu.exp.lar)) * 100## [1] 95.32425#95.32425%
otu.exp.lar.trim <- data.frame(ps.exp.lar.trim@otu_table)
sum(otu.exp.lar.trim == 0) / (ncol(otu.exp.lar.trim) * nrow(otu.exp.lar.trim)) * 100## [1] 61.81818#61.81818%
##larvae - with all
otu.exp.all.lar <- data.frame(ps.exp.all.lar@otu_table)
sum(otu.exp.all.lar == 0) / (ncol(otu.exp.all.lar) * nrow(otu.exp.all.lar)) * 100## [1] 95.96498#95.96498%
otu.exp.all.lar.trim <- data.frame(ps.exp.all.lar.trim@otu_table)
sum(otu.exp.all.lar.trim == 0) / (ncol(otu.exp.all.lar.trim) * nrow(otu.exp.all.lar.trim)) * 100## [1] 70.96296#70.96296%
##save if ready
#saveRDS(ps.exp.lar.trim,file="ps.exp.lar.trim.rds")
#saveRDS(ps.exp.all.lar.trim,file="ps.exp.all.lar.trim.rds")2.1.1 Which ASVs to subset from fasta
Didn’t do yet
otus.lar.trim <- colnames(otu.exp.lar.trim)
otus.wtr.trim <- colnames(otu.exp.wtr.trim)
#write.table(otus.lar.trim, file="otus.larvae.trim.txt", append = FALSE, sep = "/n", row.names = FALSE, col.names = FALSE,quote=FALSE)
#write.table(otus.wtr.trim, file="otus.water.trim.txt", append = FALSE, sep = "/n", row.names = FALSE, col.names = FALSE,quote=FALSE)2.2 More formatting - long form
2.2.1 Just EMO, MMO, and MEM
otu.wtr <- data.frame(ps.exp.wtr.trim@otu_table)
otu.lar <- data.frame(ps.exp.lar.trim@otu_table)
#total read count colums
otu.wtr$sum <- rowSums(otu.wtr)
otu.lar$sum <- rowSums(otu.lar)
otu.wtr$Short_label <- row.names(otu.wtr)
otu.lar$Short_label <- row.names(otu.lar)
#bring in the metadata again
samdf.wtr <- data.frame(ps.exp.wtr.trim@sam_data)
samdf.lar <- data.frame(ps.exp.lar.trim@sam_data)
##add extra metadata
#meso.data <- read.csv("../bamboomesos_metadata.csv")
#samdf.wtr.more <- merge(samdf.wtr,meso.data)
#samdf.lar.more <- merge(samdf.lar,meso.data)
#samdf.wtr$food.microbes <- paste0(samdf.wtr$Food_type,"_",samdf.wtr$Microbe_treatment)
#samdf.lar$food.microbes <- paste0(samdf.lar$Food_type,"_",samdf.lar$Microbe_treatment)
#
# #join the data sets by sample name
df.water <- plyr::join(otu.wtr, samdf.wtr, by="Short_label", type="left")
df.larv <- plyr::join(otu.lar, samdf.lar, by="Short_label", type="left")
colnames(df.water)## [1] "ASV0002" "ASV0003" "ASV0004"
## [4] "ASV0005" "ASV0007" "ASV0008"
## [7] "ASV0009" "ASV0010" "ASV0011"
## [10] "ASV0012" "ASV0015" "ASV0016"
## [13] "ASV0017" "ASV0020" "ASV0021"
## [16] "ASV0022" "ASV0023" "ASV0024"
## [19] "ASV0025" "ASV0027" "ASV0036"
## [22] "ASV0037" "ASV0039" "ASV0040"
## [25] "ASV0042" "ASV0044" "ASV0051"
## [28] "ASV0054" "ASV0057" "ASV0059"
## [31] "ASV0063" "ASV0065" "ASV0066"
## [34] "ASV0069" "ASV0072" "ASV0073"
## [37] "ASV0079" "ASV0080" "ASV0082"
## [40] "ASV0084" "ASV0085" "ASV0089"
## [43] "ASV0091" "ASV0093" "ASV0094"
## [46] "ASV0095" "ASV0101" "ASV0102"
## [49] "ASV0105" "ASV0111" "ASV0112"
## [52] "ASV0115" "ASV0116" "ASV0119"
## [55] "ASV0123" "ASV0128" "ASV0139"
## [58] "ASV0145" "ASV0146" "ASV0155"
## [61] "ASV0167" "ASV0180" "sum"
## [64] "Short_label" "Longer_name" "Sample_type"
## [67] "Num_larvae" "Date_collected" "Exp_day"
## [70] "Mesocosm_id" "Mesocosm_treatment" "Microbe_treatment"
## [73] "Food_type" "Mesocosm_type" "Treatment_notes"
## [76] "Notes" "Raw_reads" "Num_larvae_d00"
## [79] "Num_larvae_d08" "Num_larvae_d34" "Num_larvae_d40"
## [82] "Total_pupae" "Proportion_pupated" "Biomass_day30"
## [85] "Biomass_day40" "lib_size_clean" "is.neg"## [1] "ASV0002" "ASV0003" "ASV0004"
## [4] "ASV0005" "ASV0007" "ASV0008"
## [7] "ASV0009" "ASV0010" "ASV0011"
## [10] "ASV0012" "ASV0015" "ASV0016"
## [13] "ASV0017" "ASV0020" "ASV0021"
## [16] "ASV0022" "ASV0023" "ASV0024"
## [19] "ASV0025" "ASV0027" "ASV0036"
## [22] "ASV0037" "ASV0039" "ASV0040"
## [25] "ASV0042" "ASV0044" "ASV0051"
## [28] "ASV0054" "ASV0057" "ASV0059"
## [31] "ASV0063" "ASV0065" "ASV0066"
## [34] "ASV0069" "ASV0072" "ASV0073"
## [37] "ASV0079" "ASV0080" "ASV0082"
## [40] "ASV0084" "ASV0085" "ASV0089"
## [43] "ASV0091" "ASV0093" "ASV0094"
## [46] "ASV0095" "ASV0101" "ASV0102"
## [49] "ASV0105" "ASV0111" "ASV0112"
## [52] "ASV0115" "ASV0116" "ASV0119"
## [55] "ASV0123" "ASV0128" "ASV0139"
## [58] "ASV0145" "ASV0146" "ASV0155"
## [61] "ASV0167" "ASV0180" "sum"
## [64] "Short_label" "Longer_name" "Sample_type"
## [67] "Num_larvae" "Date_collected" "Exp_day"
## [70] "Mesocosm_id" "Mesocosm_treatment" "Microbe_treatment"
## [73] "Food_type" "Mesocosm_type" "Treatment_notes"
## [76] "Notes" "Raw_reads" "Num_larvae_d00"
## [79] "Num_larvae_d08" "Num_larvae_d34" "Num_larvae_d40"
## [82] "Total_pupae" "Proportion_pupated" "Biomass_day30"
## [85] "Biomass_day40" "lib_size_clean" "is.neg"df.water1 <- df.water %>%
dplyr::select(-Short_label,-Longer_name,-Sample_type,-Num_larvae,-Date_collected,-Exp_day,-Mesocosm_treatment,-Mesocosm_type,-Treatment_notes,-Notes,-Raw_reads,-Num_larvae_d00,-Num_larvae_d08,-Num_larvae_d34,-Num_larvae_d40,-Total_pupae,-Proportion_pupated,-Biomass_day30,-lib_size_clean,-is.neg)
colnames(df.water1)## [1] "ASV0002" "ASV0003" "ASV0004"
## [4] "ASV0005" "ASV0007" "ASV0008"
## [7] "ASV0009" "ASV0010" "ASV0011"
## [10] "ASV0012" "ASV0015" "ASV0016"
## [13] "ASV0017" "ASV0020" "ASV0021"
## [16] "ASV0022" "ASV0023" "ASV0024"
## [19] "ASV0025" "ASV0027" "ASV0036"
## [22] "ASV0037" "ASV0039" "ASV0040"
## [25] "ASV0042" "ASV0044" "ASV0051"
## [28] "ASV0054" "ASV0057" "ASV0059"
## [31] "ASV0063" "ASV0065" "ASV0066"
## [34] "ASV0069" "ASV0072" "ASV0073"
## [37] "ASV0079" "ASV0080" "ASV0082"
## [40] "ASV0084" "ASV0085" "ASV0089"
## [43] "ASV0091" "ASV0093" "ASV0094"
## [46] "ASV0095" "ASV0101" "ASV0102"
## [49] "ASV0105" "ASV0111" "ASV0112"
## [52] "ASV0115" "ASV0116" "ASV0119"
## [55] "ASV0123" "ASV0128" "ASV0139"
## [58] "ASV0145" "ASV0146" "ASV0155"
## [61] "ASV0167" "ASV0180" "sum"
## [64] "Mesocosm_id" "Microbe_treatment" "Food_type"
## [67] "Biomass_day40"df.water.long <- melt(df.water1, id.vars=c("sum", "Mesocosm_id", "Microbe_treatment","Food_type", "Biomass_day40"))
df.water.long$rowid <- 1:nrow(df.water.long)
df.larv1 <- df.larv %>%
dplyr::select(-Short_label,-Longer_name,-Sample_type,-Num_larvae,-Date_collected,-Exp_day,-Mesocosm_treatment,-Mesocosm_type,-Treatment_notes,-Notes,-Raw_reads,-Num_larvae_d00,-Num_larvae_d08,-Num_larvae_d34,-Num_larvae_d40,-Total_pupae,-Proportion_pupated,-Biomass_day30,-lib_size_clean,-is.neg)
colnames(df.larv1)## [1] "ASV0002" "ASV0003" "ASV0005"
## [4] "ASV0006" "ASV0007" "ASV0009"
## [7] "ASV0010" "ASV0011" "ASV0012"
## [10] "ASV0013" "ASV0016" "ASV0019"
## [13] "ASV0024" "ASV0027" "ASV0032"
## [16] "ASV0033" "ASV0043" "ASV0053"
## [19] "ASV0054" "ASV0055" "ASV0058"
## [22] "ASV0059" "ASV0060" "ASV0061"
## [25] "ASV0066" "ASV0070" "ASV0079"
## [28] "ASV0083" "ASV0085" "ASV0087"
## [31] "ASV0099" "ASV0100" "ASV0102"
## [34] "ASV0125" "ASV0130" "ASV0134"
## [37] "ASV0138" "ASV0149" "ASV0150"
## [40] "ASV0164" "ASV0218" "ASV0262"
## [43] "ASV0289" "ASV0311" "sum"
## [46] "Mesocosm_id" "Microbe_treatment" "Food_type"
## [49] "Biomass_day40"#df.larv2 <- df.larv1[df.larv1$sum!=0,]
#df.larv2$sum==0
df.larv.long <- melt(df.larv1, id.vars=c("sum", "Mesocosm_id", "Microbe_treatment","Food_type", "Biomass_day40"))
df.larv.long$rowid <- 1:nrow(df.larv.long)
##save
#saveRDS(df.water.long,file="df.water.long.rds")
#saveRDS(df.larv.long,file="df.larv.long.rds")2.2.2 All of ’em
otu.all.wtr <- data.frame(ps.exp.all.wtr.trim@otu_table)
otu.all.lar <- data.frame(ps.exp.all.lar.trim@otu_table)
#total read count colums
otu.all.wtr$sum <- rowSums(otu.all.wtr)
otu.all.lar$sum <- rowSums(otu.all.lar)
otu.all.wtr$Short_label <- row.names(otu.all.wtr)
otu.all.lar$Short_label <- row.names(otu.all.lar)
#bring in the metadata again
samdf.all.wtr <- data.frame(ps.exp.all.wtr.trim@sam_data)
samdf.all.lar <- data.frame(ps.exp.all.lar.trim@sam_data)
##add extra metadata
#meso.data <- read.csv("../bamboomesos_metadata.csv")
#samdf.wtr.more <- merge(samdf.wtr,meso.data)
#samdf.lar.more <- merge(samdf.lar,meso.data)
#samdf.wtr$food.microbes <- paste0(samdf.wtr$Food_type,"_",samdf.wtr$Microbe_treatment)
#samdf.lar$food.microbes <- paste0(samdf.lar$Food_type,"_",samdf.lar$Microbe_treatment)
#
# #join the data sets by sample name
df.all.water <- plyr::join(otu.all.wtr, samdf.all.wtr, by="Short_label", type="left")
df.all.larv <- plyr::join(otu.all.lar, samdf.all.lar, by="Short_label", type="left")
colnames(df.all.water)## [1] "ASV0001" "ASV0002" "ASV0003"
## [4] "ASV0004" "ASV0005" "ASV0007"
## [7] "ASV0008" "ASV0009" "ASV0010"
## [10] "ASV0011" "ASV0012" "ASV0015"
## [13] "ASV0016" "ASV0017" "ASV0018"
## [16] "ASV0020" "ASV0021" "ASV0022"
## [19] "ASV0023" "ASV0024" "ASV0025"
## [22] "ASV0027" "ASV0033" "ASV0036"
## [25] "ASV0037" "ASV0039" "ASV0040"
## [28] "ASV0042" "ASV0044" "ASV0051"
## [31] "ASV0054" "ASV0057" "ASV0059"
## [34] "ASV0063" "ASV0065" "ASV0066"
## [37] "ASV0069" "ASV0072" "ASV0073"
## [40] "ASV0079" "ASV0080" "ASV0082"
## [43] "ASV0083" "ASV0084" "ASV0085"
## [46] "ASV0089" "ASV0091" "ASV0093"
## [49] "ASV0094" "ASV0095" "ASV0101"
## [52] "ASV0102" "ASV0105" "ASV0111"
## [55] "ASV0112" "ASV0115" "ASV0116"
## [58] "ASV0119" "ASV0123" "ASV0128"
## [61] "ASV0139" "ASV0145" "ASV0146"
## [64] "ASV0152" "ASV0155" "ASV0167"
## [67] "ASV0175" "ASV0180" "ASV0184"
## [70] "sum" "Short_label" "Longer_name"
## [73] "Sample_type" "Num_larvae" "Date_collected"
## [76] "Exp_day" "Mesocosm_id" "Mesocosm_treatment"
## [79] "Microbe_treatment" "Food_type" "Mesocosm_type"
## [82] "Treatment_notes" "Notes" "Raw_reads"
## [85] "Num_larvae_d00" "Num_larvae_d08" "Num_larvae_d34"
## [88] "Num_larvae_d40" "Total_pupae" "Proportion_pupated"
## [91] "Biomass_day30" "Biomass_day40" "lib_size_clean"
## [94] "is.neg"##getting rid of some irrelevant metadat
df.all.water1 <- df.all.water %>%
dplyr::select(-Short_label,-Longer_name,-Sample_type,-Num_larvae,-Date_collected,-Exp_day,-Mesocosm_treatment,-Mesocosm_type,-Treatment_notes,-Notes,-Raw_reads,-Num_larvae_d00,-Num_larvae_d08,-Num_larvae_d34,-Num_larvae_d40,-Total_pupae,-Proportion_pupated,-Biomass_day30,-lib_size_clean,-is.neg)
colnames(df.all.water1)## [1] "ASV0001" "ASV0002" "ASV0003"
## [4] "ASV0004" "ASV0005" "ASV0007"
## [7] "ASV0008" "ASV0009" "ASV0010"
## [10] "ASV0011" "ASV0012" "ASV0015"
## [13] "ASV0016" "ASV0017" "ASV0018"
## [16] "ASV0020" "ASV0021" "ASV0022"
## [19] "ASV0023" "ASV0024" "ASV0025"
## [22] "ASV0027" "ASV0033" "ASV0036"
## [25] "ASV0037" "ASV0039" "ASV0040"
## [28] "ASV0042" "ASV0044" "ASV0051"
## [31] "ASV0054" "ASV0057" "ASV0059"
## [34] "ASV0063" "ASV0065" "ASV0066"
## [37] "ASV0069" "ASV0072" "ASV0073"
## [40] "ASV0079" "ASV0080" "ASV0082"
## [43] "ASV0083" "ASV0084" "ASV0085"
## [46] "ASV0089" "ASV0091" "ASV0093"
## [49] "ASV0094" "ASV0095" "ASV0101"
## [52] "ASV0102" "ASV0105" "ASV0111"
## [55] "ASV0112" "ASV0115" "ASV0116"
## [58] "ASV0119" "ASV0123" "ASV0128"
## [61] "ASV0139" "ASV0145" "ASV0146"
## [64] "ASV0152" "ASV0155" "ASV0167"
## [67] "ASV0175" "ASV0180" "ASV0184"
## [70] "sum" "Mesocosm_id" "Microbe_treatment"
## [73] "Food_type" "Biomass_day40"df.all.water.long <- melt(df.all.water1, id.vars=c("sum", "Mesocosm_id", "Microbe_treatment","Food_type", "Biomass_day40"))
df.all.water.long$rowid <- 1:nrow(df.all.water.long)
df.all.larv1 <- df.all.larv %>%
dplyr::select(-Short_label,-Longer_name,-Sample_type,-Num_larvae,-Date_collected,-Exp_day,-Mesocosm_treatment,-Mesocosm_type,-Treatment_notes,-Notes,-Raw_reads,-Num_larvae_d00,-Num_larvae_d08,-Num_larvae_d34,-Num_larvae_d40,-Total_pupae,-Proportion_pupated,-Biomass_day30,-lib_size_clean,-is.neg)
colnames(df.all.larv1)## [1] "ASV0001" "ASV0002" "ASV0003"
## [4] "ASV0004" "ASV0005" "ASV0006"
## [7] "ASV0007" "ASV0009" "ASV0010"
## [10] "ASV0011" "ASV0012" "ASV0013"
## [13] "ASV0016" "ASV0019" "ASV0024"
## [16] "ASV0027" "ASV0030" "ASV0032"
## [19] "ASV0033" "ASV0043" "ASV0053"
## [22] "ASV0054" "ASV0055" "ASV0058"
## [25] "ASV0059" "ASV0060" "ASV0061"
## [28] "ASV0066" "ASV0070" "ASV0079"
## [31] "ASV0083" "ASV0085" "ASV0087"
## [34] "ASV0099" "ASV0100" "ASV0102"
## [37] "ASV0109" "ASV0125" "ASV0130"
## [40] "ASV0134" "ASV0138" "ASV0149"
## [43] "ASV0150" "ASV0152" "ASV0164"
## [46] "ASV0210" "ASV0212" "ASV0218"
## [49] "ASV0220" "ASV0225" "ASV0262"
## [52] "ASV0276" "ASV0289" "ASV0311"
## [55] "sum" "Mesocosm_id" "Microbe_treatment"
## [58] "Food_type" "Biomass_day40"#df.all.larv2 <- df.all.larv1[df.all.larv1$sum!=0,]
#df.all.larv2$sum==0
df.all.larv.long <- melt(df.all.larv1, id.vars=c("sum", "Mesocosm_id", "Microbe_treatment","Food_type", "Biomass_day40"))
df.all.larv.long$rowid <- 1:nrow(df.all.larv.long)
##save
#saveRDS(df.all.water.long,file="df.all.water.long.rds")
#saveRDS(df.all.larv.long,file="df.all.larv.long.rds")3 Running the statz
Can start from here
setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Bamboo_mesos/Bamboo_mesos/04.mixed_models")
df.water.long <- readRDS("df.water.long.rds")
df.larv.long <- readRDS("df.larv.long.rds")
df.all.water.long <- readRDS("df.all.water.long.rds")
df.all.larv.long <- readRDS("df.all.larv.long.rds")3.1 Water
## 'data.frame': 1860 obs. of 8 variables:
## $ sum : num 101694 104508 144484 122042 88536 ...
## $ Mesocosm_id : chr "B_MMO.4" "B_MMO.5" "B_EMO.1" "B_EMO.2" ...
## $ Microbe_treatment: chr "MMO" "MMO" "EMO" "EMO" ...
## $ Food_type : chr "Bamboo" "Bamboo" "Bamboo" "Bamboo" ...
## $ Biomass_day40 : num 2.612 0.468 3.949 0 0.621 ...
## $ variable : Factor w/ 62 levels "ASV0002","ASV0003",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 26733 27492 5785 9175 5203 ...
## $ rowid : int 1 2 3 4 5 6 7 8 9 10 ...df.water.long$rowid <- as.factor(df.water.long$rowid)
df.water.long$Microbe_treatment <- as.factor(df.water.long$Microbe_treatment)
df.water.long$Food_type <- as.factor(df.water.long$Food_type)
#df.water.long$food.microbes <- as.factor(df.water.long$food.microbes)
df.water.long$Mesocosm_id <- as.factor(df.water.long$Mesocosm_id)
str(df.water.long)## 'data.frame': 1860 obs. of 8 variables:
## $ sum : num 101694 104508 144484 122042 88536 ...
## $ Mesocosm_id : Factor w/ 30 levels "B_EMO.1","B_EMO.2",..: 14 15 1 2 3 4 5 6 7 8 ...
## $ Microbe_treatment: Factor w/ 3 levels "EMO","MEM","MMO": 3 3 1 1 1 1 1 2 2 2 ...
## $ Food_type : Factor w/ 2 levels "Bamboo","Larval food": 1 1 1 1 1 1 1 1 1 1 ...
## $ Biomass_day40 : num 2.612 0.468 3.949 0 0.621 ...
## $ variable : Factor w/ 62 levels "ASV0002","ASV0003",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 26733 27492 5785 9175 5203 ...
## $ rowid : Factor w/ 1860 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...##full model
water.mod.nb1 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~variable+Microbe_treatment+Food_type,data=df.water.long)
summary(water.mod.nb1)## Family: nbinom1 ( log )
## Formula:
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type)
## Dispersion: ~variable + Microbe_treatment + Food_type
## Data: df.water.long
##
## AIC BIC logLik deviance df.resid
## 15413.4 15828.1 -7631.7 15263.4 1785
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 1.377e-06 0.001173
## variable:Microbe_treatment (Intercept) 5.699e+00 2.387351
## variable:Food_type (Intercept) 7.351e+00 2.711233
## variable:Microbe_treatment:Food_type (Intercept) 1.430e-01 0.378103
## Number of obs: 1860, groups:
## variable, 62; variable:Microbe_treatment, 186; variable:Food_type, 124; variable:Microbe_treatment:Food_type, 372
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.0382 0.4739 -12.743 < 2e-16
## Microbe_treatmentMEM -0.1270 0.4487 -0.283 0.777222
## Microbe_treatmentMMO -2.7278 0.6008 -4.540 5.62e-06
## Food_typeLarval food -2.0534 0.5393 -3.807 0.000141
## Microbe_treatmentMEM:Food_typeLarval food 0.3195 0.1856 1.722 0.085124
## Microbe_treatmentMMO:Food_typeLarval food -0.2732 0.3673 -0.744 0.456962
##
## (Intercept) ***
## Microbe_treatmentMEM
## Microbe_treatmentMMO ***
## Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food .
## Microbe_treatmentMMO:Food_typeLarval food
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.2564665 0.2975009 21.030 < 2e-16 ***
## variableASV0003 -1.5088503 0.5166781 -2.920 0.003497 **
## variableASV0004 1.8987349 0.4512348 4.208 2.58e-05 ***
## variableASV0005 1.8317243 0.4502982 4.068 4.75e-05 ***
## variableASV0007 0.5783152 0.4281055 1.351 0.176737
## variableASV0008 0.2901073 0.4510498 0.643 0.520106
## variableASV0009 -0.1788562 0.4028234 -0.444 0.657038
## variableASV0010 0.7611656 0.4562629 1.668 0.095264 .
## variableASV0011 1.3035506 0.4568341 2.853 0.004325 **
## variableASV0012 1.2093853 0.4284904 2.822 0.004766 **
## variableASV0015 0.8622593 0.4869641 1.771 0.076613 .
## variableASV0016 -0.2424918 0.4650040 -0.521 0.602030
## variableASV0017 -0.5906261 0.5638137 -1.048 0.294843
## variableASV0020 1.3479492 0.4529875 2.976 0.002923 **
## variableASV0021 0.1122259 0.7669599 0.146 0.883664
## variableASV0022 0.5214020 0.4492747 1.161 0.245828
## variableASV0023 0.6708045 0.6555633 1.023 0.306190
## variableASV0024 1.2393792 0.4658077 2.661 0.007798 **
## variableASV0025 -0.3619200 0.6020512 -0.601 0.547743
## variableASV0027 -1.0459947 0.5842380 -1.790 0.073397 .
## variableASV0036 3.6861632 0.7364961 5.005 5.59e-07 ***
## variableASV0037 2.6200524 0.9447735 2.773 0.005551 **
## variableASV0039 0.0883535 0.5842289 0.151 0.879793
## variableASV0040 -0.5393587 0.6602007 -0.817 0.413950
## variableASV0042 1.8028771 0.5040531 3.577 0.000348 ***
## variableASV0044 0.8508108 0.5877397 1.448 0.147730
## variableASV0051 -1.1642234 0.4504576 -2.585 0.009751 **
## variableASV0054 -0.3722326 0.6419720 -0.580 0.562031
## variableASV0057 -0.6032034 0.5873371 -1.027 0.304414
## variableASV0059 -0.2193775 0.5489210 -0.400 0.689413
## variableASV0063 -0.4682974 0.5026902 -0.932 0.351552
## variableASV0065 -0.4523079 0.6357944 -0.711 0.476833
## variableASV0066 0.9424386 0.7085647 1.330 0.183496
## variableASV0069 -1.9006165 0.6244191 -3.044 0.002336 **
## variableASV0072 -1.7610774 0.6163257 -2.857 0.004272 **
## variableASV0073 1.6480567 0.7696031 2.141 0.032239 *
## variableASV0079 -1.3420093 0.7566616 -1.774 0.076131 .
## variableASV0080 1.0856345 0.6405383 1.695 0.090099 .
## variableASV0082 0.4498012 0.5686437 0.791 0.428940
## variableASV0084 -0.2608202 0.5945235 -0.439 0.660876
## variableASV0085 -1.1383417 0.4505072 -2.527 0.011511 *
## variableASV0089 0.0001207 0.6510603 0.000 0.999852
## variableASV0091 0.2946653 0.7864269 0.375 0.707892
## variableASV0093 -1.8915201 0.6102519 -3.100 0.001938 **
## variableASV0094 -0.5932495 0.6157738 -0.963 0.335336
## variableASV0095 0.4885778 0.6624257 0.738 0.460783
## variableASV0101 1.2566686 0.8443749 1.488 0.136676
## variableASV0102 -0.3685071 0.6095838 -0.605 0.545496
## variableASV0105 0.8932714 0.9708123 0.920 0.357506
## variableASV0111 -1.4285033 0.5539093 -2.579 0.009910 **
## variableASV0112 0.7922081 0.7841482 1.010 0.312362
## variableASV0115 0.0603164 0.6428722 0.094 0.925250
## variableASV0116 -0.1551065 0.5273773 -0.294 0.768674
## variableASV0119 -0.2673950 0.6575854 -0.407 0.684279
## variableASV0123 -0.6364850 0.6579777 -0.967 0.333376
## variableASV0128 0.7760974 0.9009622 0.861 0.389013
## variableASV0139 1.1327518 0.8767022 1.292 0.196336
## variableASV0145 -2.6420271 0.8629926 -3.061 0.002203 **
## variableASV0146 0.0246931 0.7256392 0.034 0.972854
## variableASV0155 -1.9592978 0.6293727 -3.113 0.001851 **
## variableASV0167 -0.5380176 0.7807569 -0.689 0.490762
## variableASV0180 -3.2500265 0.6435964 -5.050 4.42e-07 ***
## Microbe_treatmentMEM 0.0107150 0.1390419 0.077 0.938573
## Microbe_treatmentMMO 2.0624551 0.2556546 8.067 7.18e-16 ***
## Food_typeLarval food -0.1704443 0.1641874 -1.038 0.299220
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1water.mod.nb2 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom2,dispformula=~variable+Microbe_treatment+Food_type,data=df.water.long)
summary(water.mod.nb2)## Family: nbinom2 ( log )
## Formula:
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type)
## Dispersion: ~variable + Microbe_treatment + Food_type
## Data: df.water.long
##
## AIC BIC logLik deviance df.resid
## 15930.5 16345.2 -7890.3 15780.5 1785
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 9.984e-07 0.0009992
## variable:Microbe_treatment (Intercept) 2.456e+01 4.9561702
## variable:Food_type (Intercept) 2.132e+01 4.6173294
## variable:Microbe_treatment:Food_type (Intercept) 1.057e-01 0.3250761
## Number of obs: 1860, groups:
## variable, 62; variable:Microbe_treatment, 186; variable:Food_type, 124; variable:Microbe_treatment:Food_type, 372
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.35573 0.87820 -7.237 4.58e-13
## Microbe_treatmentMEM -0.05179 0.91574 -0.057 0.9549
## Microbe_treatmentMMO -8.93860 1.69505 -5.273 1.34e-07
## Food_typeLarval food -3.67418 0.90354 -4.066 4.77e-05
## Microbe_treatmentMEM:Food_typeLarval food 0.30694 0.20736 1.480 0.1388
## Microbe_treatmentMMO:Food_typeLarval food -0.74827 0.43617 -1.716 0.0862
##
## (Intercept) ***
## Microbe_treatmentMEM
## Microbe_treatmentMMO ***
## Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food
## Microbe_treatmentMMO:Food_typeLarval food .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.44408 0.28413 8.602 < 2e-16 ***
## variableASV0003 1.01324 0.46818 2.164 0.030447 *
## variableASV0004 -0.83267 0.43889 -1.897 0.057801 .
## variableASV0005 -0.59799 0.43798 -1.365 0.172146
## variableASV0007 -1.29878 0.39888 -3.256 0.001130 **
## variableASV0008 -2.43096 0.43906 -5.537 3.08e-08 ***
## variableASV0009 -0.17048 0.38863 -0.439 0.660899
## variableASV0010 -0.53041 0.44308 -1.197 0.231267
## variableASV0011 -1.55520 0.42533 -3.656 0.000256 ***
## variableASV0012 -2.25799 0.40576 -5.565 2.62e-08 ***
## variableASV0015 -2.91623 0.41250 -7.070 1.55e-12 ***
## variableASV0016 0.63873 0.45136 1.415 0.157030
## variableASV0017 -2.22143 0.52617 -4.222 2.42e-05 ***
## variableASV0020 -1.77903 0.41558 -4.281 1.86e-05 ***
## variableASV0021 -3.33457 0.55482 -6.010 1.85e-09 ***
## variableASV0022 -1.72471 0.41958 -4.111 3.95e-05 ***
## variableASV0023 -2.88667 0.54278 -5.318 1.05e-07 ***
## variableASV0024 -2.23507 0.42275 -5.287 1.24e-07 ***
## variableASV0025 -0.65671 0.55776 -1.177 0.239035
## variableASV0027 -2.19496 0.55443 -3.959 7.53e-05 ***
## variableASV0036 -4.71917 0.44817 -10.530 < 2e-16 ***
## variableASV0037 -5.78499 0.52382 -11.044 < 2e-16 ***
## variableASV0039 -3.91863 0.49526 -7.912 2.53e-15 ***
## variableASV0040 -2.58160 0.55868 -4.621 3.82e-06 ***
## variableASV0042 -3.02878 0.41361 -7.323 2.43e-13 ***
## variableASV0044 -1.55150 0.53856 -2.881 0.003966 **
## variableASV0051 -2.43832 0.43082 -5.660 1.52e-08 ***
## variableASV0054 -3.98059 0.45992 -8.655 < 2e-16 ***
## variableASV0057 -3.48125 0.48176 -7.226 4.97e-13 ***
## variableASV0059 -3.74888 0.45429 -8.252 < 2e-16 ***
## variableASV0063 -2.99035 0.44526 -6.716 1.87e-11 ***
## variableASV0065 -3.32762 0.51256 -6.492 8.46e-11 ***
## variableASV0066 -4.30290 0.48849 -8.809 < 2e-16 ***
## variableASV0069 0.27072 0.57313 0.472 0.636674
## variableASV0072 -0.13073 0.56912 -0.230 0.818321
## variableASV0073 -4.73197 0.49269 -9.604 < 2e-16 ***
## variableASV0079 -2.66730 0.57854 -4.610 4.02e-06 ***
## variableASV0080 -4.13395 0.45285 -9.129 < 2e-16 ***
## variableASV0082 -3.69105 0.43927 -8.403 < 2e-16 ***
## variableASV0084 -3.01801 0.47206 -6.393 1.62e-10 ***
## variableASV0085 -1.19092 0.42990 -2.770 0.005602 **
## variableASV0089 -2.62620 0.55567 -4.726 2.29e-06 ***
## variableASV0091 -3.87883 0.57103 -6.793 1.10e-11 ***
## variableASV0093 -0.61818 0.56392 -1.096 0.272979
## variableASV0094 -0.64417 0.56710 -1.136 0.256000
## variableASV0095 -2.91679 0.53455 -5.457 4.85e-08 ***
## variableASV0101 -4.77505 0.49073 -9.731 < 2e-16 ***
## variableASV0102 -2.95578 0.48183 -6.135 8.54e-10 ***
## variableASV0105 -4.40262 0.61297 -7.182 6.85e-13 ***
## variableASV0111 -3.66298 0.46197 -7.929 2.21e-15 ***
## variableASV0112 -3.52209 0.54224 -6.495 8.28e-11 ***
## variableASV0115 -4.23469 0.49012 -8.640 < 2e-16 ***
## variableASV0116 -3.54838 0.42923 -8.267 < 2e-16 ***
## variableASV0119 -2.56979 0.56030 -4.586 4.51e-06 ***
## variableASV0123 -2.68720 0.55503 -4.842 1.29e-06 ***
## variableASV0128 -3.90753 0.56794 -6.880 5.98e-12 ***
## variableASV0139 -3.97281 0.56553 -7.025 2.14e-12 ***
## variableASV0145 -2.14679 0.61828 -3.472 0.000516 ***
## variableASV0146 -4.27875 0.48611 -8.802 < 2e-16 ***
## variableASV0155 -0.09976 0.57780 -0.173 0.862922
## variableASV0167 -3.60614 0.58520 -6.162 7.17e-10 ***
## variableASV0180 0.16370 0.59327 0.276 0.782605
## Microbe_treatmentMEM -0.18359 0.12110 -1.516 0.129514
## Microbe_treatmentMMO -1.38803 0.25430 -5.458 4.81e-08 ***
## Food_typeLarval food -0.20719 0.14201 -1.459 0.144567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## water.mod.nb1 0.0 75
## water.mod.nb2 517.1 75
##dispersal formula checks
water.mod.nb1.disnovar <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~Microbe_treatment+Food_type,data=df.water.long)
anova(water.mod.nb1.disnovar,water.mod.nb1) #sig***## Data: df.water.long
## Models:
## water.mod.nb1.disnovar: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1.disnovar: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.disnovar: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## water.mod.nb1.disnovar 14 15642 15720 -7807.1 15614
## water.mod.nb1 75 15413 15828 -7631.7 15263 350.86 61 < 2.2e-16
##
## water.mod.nb1.disnovar
## water.mod.nb1 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##removing 3 way interaction
water.mod.nb1.nofm <- update(water.mod.nb1,.~.-(1|variable:Microbe_treatment:Food_type))
anova(water.mod.nb1,water.mod.nb1.nofm)## Data: df.water.long
## Models:
## water.mod.nb1.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.nofm: variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.nofm: Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## water.mod.nb1.nofm 74 15425 15834 -7638.4 15277
## water.mod.nb1 75 15413 15828 -7631.7 15263 13.473 1 0.000242 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##super sig 0.000242
water.mod.nb2.nofm <- update(water.mod.nb2,.~.-(1|variable:Microbe_treatment:Food_type))
anova(water.mod.nb2,water.mod.nb2.nofm)## Data: df.water.long
## Models:
## water.mod.nb2.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2.nofm: variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2.nofm: Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## water.mod.nb2.nofm 74 15937 16346 -7894.4 15789
## water.mod.nb2 75 15930 16345 -7890.3 15780 8.3132 1 0.003936 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##also super sig
###save conditional modes
water.mod.ranef <- as.data.frame(ranef(water.mod.nb1))
#write.csv(water.mod.ranef, "water.mod.ranef.csv")
#saveRDS(water.mod.nb1,file="water.mod.nb1.rds")3.1.1 Nbinom1 vs. nbinom2 things
Following this wonderful resource, on page 16
Ben Bolker, Mollie Brooks, Beth Gardner, Cleridy Lennert, Mihoko Minami, October 23, 2012, Owls example: a zero-inflated, generalized linear mixed model for count data. https://groups.nceas.ucsb.edu/non-linear-modeling/projects/owls/WRITEUP/owls.pdf
mvtab <- ddply(df.water.long,
.(variable:Microbe_treatment:Food_type),
summarise,
callmean=mean(value),
callvar=var(value))
q1 <- qplot(callmean,callvar,data=mvtab)## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(q1+
## linear (quasi-Poisson/NB1) fit
geom_smooth(method="lm",formula=y~x-1,colour="pink")+
## smooth (loess)
geom_smooth(colour="red")+
## semi-quadratic (NB2/LNP)
geom_smooth(method="lm",formula=y~I(x^2)+offset(x)-1,colour="purple")+
## Poisson (v=m)
geom_abline(a=0,b=1,lty=2))## Warning in geom_abline(a = 0, b = 1, lty = 2): Ignoring unknown parameters: `a`
## and `b`## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
3.2 Larvae
## 'data.frame': 1320 obs. of 8 variables:
## $ sum : num 63543 18211 13437 48065 39128 ...
## $ Mesocosm_id : chr "B_EMO.3" "B_EMO.4" "B_EMO.5" "B_MEM.1" ...
## $ Microbe_treatment: chr "EMO" "EMO" "EMO" "MEM" ...
## $ Food_type : chr "Bamboo" "Bamboo" "Bamboo" "Bamboo" ...
## $ Biomass_day40 : num 0.621 0.308 4.085 0.227 4.007 ...
## $ variable : Factor w/ 44 levels "ASV0002","ASV0003",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 0 20 0 0 0 42 0 0 257 0 ...
## $ rowid : int 1 2 3 4 5 6 7 8 9 10 ...df.larv.long$rowid <- as.factor(df.larv.long$rowid)
df.larv.long$Microbe_treatment <- as.factor(df.larv.long$Microbe_treatment)
df.larv.long$Food_type <- as.factor(df.larv.long$Food_type)
#df.larv.long$food.microbes <- as.factor(df.larv.long$food.microbes)
df.larv.long$Mesocosm_id <- as.factor(df.larv.long$Mesocosm_id)
str(df.larv.long)## 'data.frame': 1320 obs. of 8 variables:
## $ sum : num 63543 18211 13437 48065 39128 ...
## $ Mesocosm_id : Factor w/ 30 levels "B_EMO.1","B_EMO.2",..: 3 4 5 6 7 8 9 10 11 12 ...
## $ Microbe_treatment: Factor w/ 3 levels "EMO","MEM","MMO": 1 1 1 2 2 2 2 2 3 3 ...
## $ Food_type : Factor w/ 2 levels "Bamboo","Larval food": 1 1 1 1 1 1 1 1 1 1 ...
## $ Biomass_day40 : num 0.621 0.308 4.085 0.227 4.007 ...
## $ variable : Factor w/ 44 levels "ASV0002","ASV0003",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 0 20 0 0 0 42 0 0 257 0 ...
## $ rowid : Factor w/ 1320 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...##binomial
larv.mod.bin <- glmmTMB(cbind(value, sum-value)~Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type)+(1|rowid),family=binomial,data=df.larv.long)
summary(larv.mod.bin)## Family: binomial ( logit )
## Formula:
## cbind(value, sum - value) ~ Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type) + (1 | rowid)
## Data: df.larv.long
##
## AIC BIC logLik deviance df.resid
## 9953.6 10010.7 -4965.8 9931.6 1309
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 6.837e-07 0.0008269
## variable:Microbe_treatment (Intercept) 2.057e+01 4.5355850
## variable:Food_type (Intercept) 3.687e+01 6.0717861
## variable:Microbe_treatment:Food_type (Intercept) 1.325e+01 3.6401214
## rowid (Intercept) 1.003e+01 3.1666214
## Number of obs: 1320, groups:
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88; variable:Microbe_treatment:Food_type, 264; rowid, 1320
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) -14.040307 1.440472 -9.747
## Microbe_treatmentMEM 3.377608 1.410992 2.394
## Microbe_treatmentMMO -4.446516 1.613997 -2.755
## Food_typeLarval food 0.003874 1.729333 0.002
## Microbe_treatmentMEM:Food_typeLarval food -3.474420 1.498526 -2.319
## Microbe_treatmentMMO:Food_typeLarval food -2.067792 1.602383 -1.290
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## Microbe_treatmentMEM 0.01668 *
## Microbe_treatmentMMO 0.00587 **
## Food_typeLarval food 0.99821
## Microbe_treatmentMEM:Food_typeLarval food 0.02042 *
## Microbe_treatmentMMO:Food_typeLarval food 0.19689
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##full model
larv.mod.nb1 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),disp=~variable+Food_type+Microbe_treatment,family=nbinom1,data=df.larv.long)
summary(larv.mod.nb1)## Family: nbinom1 ( log )
## Formula:
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type)
## Dispersion: ~variable + Food_type + Microbe_treatment
## Data: df.larv.long
##
## AIC BIC logLik deviance df.resid
## 9374.9 9670.5 -4630.5 9260.9 1263
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 1.856e-07 0.0004308
## variable:Microbe_treatment (Intercept) 3.439e+00 1.8544060
## variable:Food_type (Intercept) 7.217e+00 2.6865139
## variable:Microbe_treatment:Food_type (Intercept) 6.505e-01 0.8065311
## Number of obs: 1320, groups:
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88; variable:Microbe_treatment:Food_type, 264
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.95183 0.58047 -11.976 < 2e-16
## Microbe_treatmentMEM 0.29763 0.50828 0.586 0.55817
## Microbe_treatmentMMO -1.62053 0.59085 -2.743 0.00609
## Food_typeLarval food 0.17975 0.70299 0.256 0.79819
## Microbe_treatmentMEM:Food_typeLarval food -0.53532 0.42755 -1.252 0.21055
## Microbe_treatmentMMO:Food_typeLarval food -0.06406 0.52157 -0.123 0.90225
##
## (Intercept) ***
## Microbe_treatmentMEM
## Microbe_treatmentMMO **
## Food_typeLarval food
## Microbe_treatmentMEM:Food_typeLarval food
## Microbe_treatmentMMO:Food_typeLarval food
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.56954 0.58181 9.573 < 2e-16 ***
## variableASV0003 2.42329 0.62861 3.855 0.000116 ***
## variableASV0005 1.09329 0.85231 1.283 0.199582
## variableASV0006 2.70004 0.65837 4.101 4.11e-05 ***
## variableASV0007 1.55254 0.66242 2.344 0.019092 *
## variableASV0009 1.99002 0.69265 2.873 0.004065 **
## variableASV0010 0.91182 0.73202 1.246 0.212905
## variableASV0011 -0.13271 0.75148 -0.177 0.859820
## variableASV0012 1.80874 1.01104 1.789 0.073615 .
## variableASV0013 2.91827 0.76372 3.821 0.000133 ***
## variableASV0016 1.53602 0.85862 1.789 0.073625 .
## variableASV0019 2.82165 0.71026 3.973 7.11e-05 ***
## variableASV0024 1.29458 0.70167 1.845 0.065037 .
## variableASV0027 1.97844 0.72858 2.715 0.006618 **
## variableASV0032 2.00643 0.79047 2.538 0.011140 *
## variableASV0033 2.64198 0.71766 3.681 0.000232 ***
## variableASV0043 4.43006 0.87997 5.034 4.79e-07 ***
## variableASV0053 1.16523 0.72785 1.601 0.109397
## variableASV0054 1.68061 0.89688 1.874 0.060952 .
## variableASV0055 2.38325 0.78221 3.047 0.002313 **
## variableASV0058 2.93112 0.84328 3.476 0.000509 ***
## variableASV0059 0.60771 0.70774 0.859 0.390530
## variableASV0060 3.10043 0.83244 3.725 0.000196 ***
## variableASV0061 1.71121 0.81456 2.101 0.035660 *
## variableASV0066 -0.47107 0.73958 -0.637 0.524162
## variableASV0070 1.37715 0.92483 1.489 0.136467
## variableASV0079 0.76562 0.68258 1.122 0.262005
## variableASV0083 0.62041 0.68273 0.909 0.363494
## variableASV0085 1.21442 1.03035 1.179 0.238536
## variableASV0087 1.81063 0.76428 2.369 0.017833 *
## variableASV0099 1.67867 0.82274 2.040 0.041317 *
## variableASV0100 1.00298 0.74743 1.342 0.179630
## variableASV0102 -0.43747 0.85386 -0.512 0.608412
## variableASV0125 2.27182 0.87847 2.586 0.009707 **
## variableASV0130 0.98610 0.71145 1.386 0.165730
## variableASV0134 2.22981 1.00464 2.220 0.026452 *
## variableASV0138 0.41231 0.88483 0.466 0.641228
## variableASV0149 -1.11486 0.70938 -1.572 0.116043
## variableASV0150 0.06528 0.81195 0.080 0.935924
## variableASV0164 -0.18252 0.73054 -0.250 0.802708
## variableASV0218 0.10549 0.94527 0.112 0.911145
## variableASV0262 0.74410 1.03491 0.719 0.472138
## variableASV0289 -0.15667 0.86232 -0.182 0.855832
## variableASV0311 0.23946 0.88198 0.272 0.786003
## Food_typeLarval food 0.93604 0.28051 3.337 0.000847 ***
## Microbe_treatmentMEM -0.77455 0.23725 -3.265 0.001096 **
## Microbe_treatmentMMO 0.05033 0.31572 0.159 0.873344
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#saveRDS(larv.mod.nb1,file="larv.mod.nb1.rds")
larv.mod.nb2 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),disp=~variable+Food_type+Microbe_treatment,family=nbinom2,data=df.larv.long)
summary(larv.mod.nb2)## Family: nbinom2 ( log )
## Formula:
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type)
## Dispersion: ~variable + Food_type + Microbe_treatment
## Data: df.larv.long
##
## AIC BIC logLik deviance df.resid
## 9752.6 10048.2 -4819.3 9638.6 1263
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 7.414e-07 0.000861
## variable:Microbe_treatment (Intercept) 1.865e+01 4.318578
## variable:Food_type (Intercept) 2.931e+01 5.413738
## variable:Microbe_treatment:Food_type (Intercept) 3.068e+00 1.751673
## Number of obs: 1320, groups:
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88; variable:Microbe_treatment:Food_type, 264
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.53380 1.18259 -8.062 7.52e-16
## Microbe_treatmentMEM 2.60767 1.11997 2.328 0.019894
## Microbe_treatmentMMO -2.67102 1.25956 -2.121 0.033955
## Food_typeLarval food 0.01895 1.37397 0.014 0.988998
## Microbe_treatmentMEM:Food_typeLarval food -2.99036 0.89874 -3.327 0.000877
## Microbe_treatmentMMO:Food_typeLarval food -2.69280 1.10875 -2.429 0.015154
##
## (Intercept) ***
## Microbe_treatmentMEM *
## Microbe_treatmentMMO *
## Food_typeLarval food
## Microbe_treatmentMEM:Food_typeLarval food ***
## Microbe_treatmentMMO:Food_typeLarval food *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.779935 0.372728 -7.458 8.76e-14 ***
## variableASV0003 1.636689 0.442891 3.695 0.000219 ***
## variableASV0005 0.518260 0.591254 0.877 0.380735
## variableASV0006 1.595294 0.496609 3.212 0.001316 **
## variableASV0007 1.145117 0.485977 2.356 0.018457 *
## variableASV0009 0.830116 0.478662 1.734 0.082875 .
## variableASV0010 0.653609 0.485088 1.347 0.177851
## variableASV0011 2.628128 0.585580 4.488 7.19e-06 ***
## variableASV0012 -0.861885 0.648264 -1.330 0.183674
## variableASV0013 4.181151 0.598683 6.984 2.87e-12 ***
## variableASV0016 0.487424 0.590490 0.825 0.409112
## variableASV0019 1.813278 0.519091 3.493 0.000477 ***
## variableASV0024 0.947252 0.494014 1.917 0.055179 .
## variableASV0027 1.446153 0.526480 2.747 0.006017 **
## variableASV0032 3.822804 0.592276 6.454 1.09e-10 ***
## variableASV0033 0.702789 0.462448 1.520 0.128583
## variableASV0043 0.743267 0.610073 1.218 0.223100
## variableASV0053 0.945001 0.503333 1.877 0.060451 .
## variableASV0054 -0.804134 0.527172 -1.525 0.127166
## variableASV0055 2.566564 0.558195 4.598 4.27e-06 ***
## variableASV0058 1.378698 0.593636 2.322 0.020208 *
## variableASV0059 0.975630 0.512528 1.904 0.056967 .
## variableASV0060 0.513368 0.536041 0.958 0.338213
## variableASV0061 0.474951 0.499538 0.951 0.341716
## variableASV0066 1.136746 0.593298 1.916 0.055368 .
## variableASV0070 1.898357 0.628794 3.019 0.002536 **
## variableASV0079 1.510157 0.478858 3.154 0.001612 **
## variableASV0083 1.551743 0.480928 3.227 0.001253 **
## variableASV0085 -0.545332 0.610078 -0.894 0.371390
## variableASV0087 2.811920 0.564989 4.977 6.46e-07 ***
## variableASV0099 -0.396072 0.495977 -0.799 0.424541
## variableASV0100 0.979286 0.547645 1.788 0.073747 .
## variableASV0102 -0.494742 0.528771 -0.936 0.349455
## variableASV0125 0.788877 0.547143 1.442 0.149356
## variableASV0130 0.169561 0.454711 0.373 0.709224
## variableASV0134 -0.822009 0.552742 -1.487 0.136976
## variableASV0138 1.647222 0.599422 2.748 0.005996 **
## variableASV0149 3.139577 0.557746 5.629 1.81e-08 ***
## variableASV0150 1.208037 0.582530 2.074 0.038100 *
## variableASV0164 0.800755 0.511426 1.566 0.117412
## variableASV0218 0.323392 0.630189 0.513 0.607835
## variableASV0262 0.118794 0.672205 0.177 0.859726
## variableASV0289 1.829776 0.615928 2.971 0.002971 **
## variableASV0311 1.711515 0.605002 2.829 0.004670 **
## Food_typeLarval food 0.884808 0.167399 5.286 1.25e-07 ***
## Microbe_treatmentMEM 0.251416 0.168605 1.491 0.135922
## Microbe_treatmentMMO -0.008739 0.223269 -0.039 0.968779
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## larv.mod.nb1 0.0 57
## larv.mod.nb2 377.7 57
## larv.mod.bin 578.7 11## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##dispersal formula checks
larv.mod.nb1.disnovar <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~Microbe_treatment+Food_type,data=df.larv.long)
anova(larv.mod.nb1.disnovar,larv.mod.nb1) #sig***## Data: df.larv.long
## Models:
## larv.mod.nb1.disnovar: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1.disnovar: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.disnovar: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
## Df AIC BIC logLik deviance Chisq Chi Df
## larv.mod.nb1.disnovar 14 9456.8 9529.4 -4714.4 9428.8
## larv.mod.nb1 57 9374.9 9670.5 -4630.5 9260.9 167.92 43
## Pr(>Chisq)
## larv.mod.nb1.disnovar
## larv.mod.nb1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##removing 3 way interaction
larv.mod.nb1.nofm <- update(larv.mod.nb1,.~.-(1|variable:Microbe_treatment:Food_type))
summary(larv.mod.nb1.nofm)## Family: nbinom1 ( log )
## Formula:
## value ~ Microbe_treatment + Food_type + (1 | variable) + (1 |
## variable:Microbe_treatment) + (1 | variable:Food_type) +
## Microbe_treatment:Food_type + offset(log(sum))
## Dispersion: ~variable + Food_type + Microbe_treatment
## Data: df.larv.long
##
## AIC BIC logLik deviance df.resid
## 9380.4 9670.8 -4634.2 9268.4 1264
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 1.833e-07 0.0004282
## variable:Microbe_treatment (Intercept) 4.033e+00 2.0082238
## variable:Food_type (Intercept) 7.592e+00 2.7553900
## Number of obs: 1320, groups:
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.88584 0.57879 -11.897 < 2e-16
## Microbe_treatmentMEM 0.34474 0.48995 0.704 0.48167
## Microbe_treatmentMMO -1.70112 0.56011 -3.037 0.00239
## Food_typeLarval food 0.05614 0.68107 0.082 0.93431
## Microbe_treatmentMEM:Food_typeLarval food -0.47751 0.28648 -1.667 0.09555
## Microbe_treatmentMMO:Food_typeLarval food 0.04079 0.33684 0.121 0.90361
##
## (Intercept) ***
## Microbe_treatmentMEM
## Microbe_treatmentMMO **
## Food_typeLarval food
## Microbe_treatmentMEM:Food_typeLarval food .
## Microbe_treatmentMMO:Food_typeLarval food
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.66150 0.58618 9.658 < 2e-16 ***
## variableASV0003 2.32708 0.63327 3.675 0.000238 ***
## variableASV0005 1.09430 0.85645 1.278 0.201352
## variableASV0006 2.73578 0.66927 4.088 4.36e-05 ***
## variableASV0007 1.53363 0.66346 2.312 0.020803 *
## variableASV0009 1.95568 0.69320 2.821 0.004784 **
## variableASV0010 0.86886 0.73524 1.182 0.237309
## variableASV0011 -0.13682 0.75431 -0.181 0.856065
## variableASV0012 1.79164 1.01061 1.773 0.076255 .
## variableASV0013 2.81004 0.76910 3.654 0.000259 ***
## variableASV0016 1.53770 0.86280 1.782 0.074714 .
## variableASV0019 2.77594 0.71214 3.898 9.70e-05 ***
## variableASV0024 1.27234 0.70603 1.802 0.071526 .
## variableASV0027 1.88812 0.73373 2.573 0.010073 *
## variableASV0032 1.93115 0.79952 2.415 0.015718 *
## variableASV0033 2.63619 0.72605 3.631 0.000282 ***
## variableASV0043 4.33203 0.88658 4.886 1.03e-06 ***
## variableASV0053 1.06864 0.73318 1.458 0.144966
## variableASV0054 1.78701 0.92632 1.929 0.053713 .
## variableASV0055 2.28938 0.78728 2.908 0.003638 **
## variableASV0058 2.83787 0.84856 3.344 0.000825 ***
## variableASV0059 0.57687 0.71151 0.811 0.417494
## variableASV0060 3.00739 0.83668 3.594 0.000325 ***
## variableASV0061 1.60711 0.81783 1.965 0.049405 *
## variableASV0066 -0.48146 0.74257 -0.648 0.516746
## variableASV0070 1.35318 0.94334 1.434 0.151441
## variableASV0079 0.66061 0.68510 0.964 0.334922
## variableASV0083 0.51392 0.68511 0.750 0.453175
## variableASV0085 1.19338 1.03223 1.156 0.247630
## variableASV0087 1.69651 0.76889 2.206 0.027353 *
## variableASV0099 1.64988 0.82800 1.993 0.046305 *
## variableASV0100 0.99953 0.75172 1.330 0.183631
## variableASV0102 -0.41224 0.86194 -0.478 0.632457
## variableASV0125 2.14935 0.88169 2.438 0.014778 *
## variableASV0130 1.17259 0.73248 1.601 0.109411
## variableASV0134 2.21104 1.01097 2.187 0.028740 *
## variableASV0138 0.30410 0.88967 0.342 0.732493
## variableASV0149 -1.22305 0.71326 -1.715 0.086397 .
## variableASV0150 0.04696 0.81330 0.058 0.953952
## variableASV0164 -0.19801 0.73680 -0.269 0.788125
## variableASV0218 0.07840 0.94312 0.083 0.933754
## variableASV0262 0.71369 1.03328 0.691 0.489751
## variableASV0289 -0.25622 0.86861 -0.295 0.768007
## variableASV0311 0.13310 0.88717 0.150 0.880744
## Food_typeLarval food 0.84689 0.27662 3.062 0.002202 **
## Microbe_treatmentMEM -0.74771 0.23252 -3.216 0.001302 **
## Microbe_treatmentMMO 0.05841 0.31114 0.188 0.851087
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.larv.long
## Models:
## larv.mod.nb1.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.nofm: variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.nofm: Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## larv.mod.nb1.nofm 56 9380.4 9670.8 -4634.2 9268.4
## larv.mod.nb1 57 9374.9 9670.5 -4630.5 9260.9 7.5229 1 0.006092 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ##removing just food
# larv.mod.nb1.nofood <- update(larv.mod.nb1,.~.-(1|variable:Food_type))
# anova(larv.mod.nb1,larv.mod.nb1.nofood) #super sig***
#
# ##removing just microbes
# larv.mod.nb1.nomicr <- update(larv.mod.nb1,.~.-(1|variable:Microbe_treatment))
# anova(larv.mod.nb1,larv.mod.nb1.nomicr) #super sig***
##checking out nbinom2
larv.mod.nb2.nofm <- update(larv.mod.nb2,.~.-(1|variable:Microbe_treatment:Food_type))
anova(larv.mod.nb2,larv.mod.nb2.nofm)## Data: df.larv.long
## Models:
## larv.mod.nb2.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2.nofm: variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2.nofm: Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## larv.mod.nb2.nofm 56 9763.9 10054 -4825.9 9651.9
## larv.mod.nb2 57 9752.6 10048 -4819.3 9638.6 13.281 1 0.0002681 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##sig**
##checking out binomial
larv.mod.bin.nofm <- update(larv.mod.bin,.~.-(1|variable:Microbe_treatment:Food_type))
summary(larv.mod.bin.nofm)## Family: binomial ( logit )
## Formula:
## cbind(value, sum - value) ~ Microbe_treatment + Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | rowid) + Microbe_treatment:Food_type
## Data: df.larv.long
##
## AIC BIC logLik deviance df.resid
## 9979.2 10031.1 -4979.6 9959.2 1310
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 9.084e-07 0.0009531
## variable:Microbe_treatment (Intercept) 3.206e+01 5.6618257
## variable:Food_type (Intercept) 4.579e+01 6.7668655
## rowid (Intercept) 1.071e+01 3.2732847
## Number of obs: 1320, groups:
## variable, 44; variable:Microbe_treatment, 132; variable:Food_type, 88; rowid, 1320
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -13.4410 1.4590 -9.213 < 2e-16
## Microbe_treatmentMEM 3.0963 1.3187 2.348 0.018878
## Microbe_treatmentMMO -5.6255 1.4865 -3.784 0.000154
## Food_typeLarval food -0.5707 1.6528 -0.345 0.729870
## Microbe_treatmentMEM:Food_typeLarval food -2.2886 0.7195 -3.181 0.001470
## Microbe_treatmentMMO:Food_typeLarval food -0.9151 0.9024 -1.014 0.310532
##
## (Intercept) ***
## Microbe_treatmentMEM *
## Microbe_treatmentMMO ***
## Food_typeLarval food
## Microbe_treatmentMEM:Food_typeLarval food **
## Microbe_treatmentMMO:Food_typeLarval food
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.larv.long
## Models:
## larv.mod.bin.nofm: cbind(value, sum - value) ~ Microbe_treatment + Food_type + (1 | , zi=~0, disp=~1
## larv.mod.bin.nofm: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~1
## larv.mod.bin.nofm: (1 | rowid) + Microbe_treatment:Food_type, zi=~0, disp=~1
## larv.mod.bin: cbind(value, sum - value) ~ Microbe_treatment * Food_type + (1 | , zi=~0, disp=~1
## larv.mod.bin: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~1
## larv.mod.bin: (1 | variable:Microbe_treatment:Food_type) + (1 | rowid), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## larv.mod.bin.nofm 10 9979.2 10031 -4979.6 9959.2
## larv.mod.bin 11 9953.6 10011 -4965.8 9931.6 27.627 1 1.471e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1###save conditional modes
larv.mod.ranef <- as.data.frame(ranef(larv.mod.nb1))
#write.csv(larv.mod.ranef, "larv.mod.ranef.csv")
# library("ggstats")
# ggcoef_multicomponents(larv.mod.nb1)3.2.1 Nbinom1 vs. nbinom2 things
##note: some errors/warnings if the treatments aren't factors for some reason
mvtab <- ddply(df.larv.long,
.(variable:Microbe_treatment:Food_type),
summarise,
callmean=mean(value),
callvar=var(value))
q1 <- qplot(callmean,callvar,data=mvtab)
print(q1+
## linear (quasi-Poisson/NB1) fit
geom_smooth(method="lm",formula=y~x-1,colour="pink")+
## smooth (loess)
geom_smooth(colour="red")+
## semi-quadratic (NB2/LNP)
geom_smooth(method="lm",formula=y~I(x^2)+offset(x)-1,colour="purple")+
## Poisson (v=m)
geom_abline(a=0,b=1,lty=2))#+## Warning in geom_abline(a = 0, b = 1, lty = 2): Ignoring unknown parameters: `a`
## and `b`## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
3.3 Water - all treatments
## 'data.frame': 3450 obs. of 8 variables:
## $ sum : num 101895 104735 144484 122710 89143 ...
## $ Mesocosm_id : chr "B_MMO.4" "B_MMO.5" "B_EMO.1" "B_EMO.2" ...
## $ Microbe_treatment: chr "MMO" "MMO" "EMO" "EMO" ...
## $ Food_type : chr "Bamboo" "Bamboo" "Bamboo" "Bamboo" ...
## $ Biomass_day40 : num 2.612 0.468 3.949 0 0.621 ...
## $ variable : Factor w/ 69 levels "ASV0001","ASV0002",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 0 0 0 0 0 ...
## $ rowid : int 1 2 3 4 5 6 7 8 9 10 ...df.all.water.long$rowid <- as.factor(df.all.water.long$rowid)
df.all.water.long$Microbe_treatment <- as.factor(df.all.water.long$Microbe_treatment)
df.all.water.long$Food_type <- as.factor(df.all.water.long$Food_type)
#df.all.water.long$food.microbes <- as.factor(df.all.water.long$food.microbes)
df.all.water.long$Mesocosm_id <- as.factor(df.all.water.long$Mesocosm_id)
str(df.all.water.long)## 'data.frame': 3450 obs. of 8 variables:
## $ sum : num 101895 104735 144484 122710 89143 ...
## $ Mesocosm_id : Factor w/ 50 levels "B_ALL.1","B_ALL.2",..: 24 25 11 12 13 14 15 6 7 8 ...
## $ Microbe_treatment: Factor w/ 5 levels "ALL","ECO","EMO",..: 5 5 3 3 3 3 3 2 2 2 ...
## $ Food_type : Factor w/ 2 levels "Bamboo","Larval food": 1 1 1 1 1 1 1 1 1 1 ...
## $ Biomass_day40 : num 2.612 0.468 3.949 0 0.621 ...
## $ variable : Factor w/ 69 levels "ASV0001","ASV0002",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 0 0 0 0 0 ...
## $ rowid : Factor w/ 3450 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...##full model
water.mod.nb1 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~variable+Microbe_treatment+Food_type,data=df.all.water.long)
summary(water.mod.nb1)## Family: nbinom1 ( log )
## Formula:
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type)
## Dispersion: ~variable + Microbe_treatment + Food_type
## Data: df.all.water.long
##
## AIC BIC logLik deviance df.resid
## 22998.2 23539.0 -11411.1 22822.2 3362
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 0.2286 0.4781
## variable:Microbe_treatment (Intercept) 6.3695 2.5238
## variable:Food_type (Intercept) 7.8079 2.7943
## variable:Microbe_treatment:Food_type (Intercept) 0.2461 0.4961
## Number of obs: 3450, groups:
## variable, 69; variable:Microbe_treatment, 345; variable:Food_type, 138; variable:Microbe_treatment:Food_type, 690
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.44631 0.47687 -13.518 < 2e-16
## Microbe_treatmentECO -7.58589 0.80260 -9.452 < 2e-16
## Microbe_treatmentEMO 0.08668 0.45898 0.189 0.850211
## Microbe_treatmentMEM -0.12823 0.45788 -0.280 0.779435
## Microbe_treatmentMMO -2.51900 0.57718 -4.364 1.28e-05
## Food_typeLarval food -2.66398 0.54486 -4.889 1.01e-06
## Microbe_treatmentECO:Food_typeLarval food -1.34724 0.58473 -2.304 0.021220
## Microbe_treatmentEMO:Food_typeLarval food 0.56380 0.23769 2.372 0.017691
## Microbe_treatmentMEM:Food_typeLarval food 0.90291 0.23485 3.845 0.000121
## Microbe_treatmentMMO:Food_typeLarval food 0.54351 0.37460 1.451 0.146814
##
## (Intercept) ***
## Microbe_treatmentECO ***
## Microbe_treatmentEMO
## Microbe_treatmentMEM
## Microbe_treatmentMMO ***
## Food_typeLarval food ***
## Microbe_treatmentECO:Food_typeLarval food *
## Microbe_treatmentEMO:Food_typeLarval food *
## Microbe_treatmentMEM:Food_typeLarval food ***
## Microbe_treatmentMMO:Food_typeLarval food
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.37016 0.42392 15.027 < 2e-16 ***
## variableASV0002 0.77344 0.46332 1.669 0.095052 .
## variableASV0003 -0.92224 0.55333 -1.667 0.095571 .
## variableASV0004 2.63337 0.51177 5.146 2.67e-07 ***
## variableASV0005 2.25959 0.50658 4.461 8.18e-06 ***
## variableASV0007 0.97636 0.50479 1.934 0.053090 .
## variableASV0008 1.31909 0.51861 2.544 0.010974 *
## variableASV0009 0.21037 0.48832 0.431 0.666613
## variableASV0010 1.43391 0.50300 2.851 0.004362 **
## variableASV0011 1.88834 0.51543 3.664 0.000249 ***
## variableASV0012 1.49855 0.50859 2.947 0.003214 **
## variableASV0015 2.47173 0.55981 4.415 1.01e-05 ***
## variableASV0016 0.08453 0.51955 0.163 0.870757
## variableASV0017 0.09542 0.58335 0.164 0.870062
## variableASV0018 -0.11272 0.89425 -0.126 0.899692
## variableASV0020 1.87992 0.53161 3.536 0.000406 ***
## variableASV0021 1.52117 0.76308 1.993 0.046210 *
## variableASV0022 1.07533 0.51278 2.097 0.035988 *
## variableASV0023 1.46196 0.61914 2.361 0.018212 *
## variableASV0024 1.40504 0.52409 2.681 0.007342 **
## variableASV0025 -0.21316 0.60895 -0.350 0.726302
## variableASV0027 -0.73196 0.59584 -1.228 0.219277
## variableASV0033 1.39685 0.96018 1.455 0.145731
## variableASV0036 3.75485 0.67649 5.550 2.85e-08 ***
## variableASV0037 3.54829 0.86838 4.086 4.39e-05 ***
## variableASV0039 0.73133 0.58984 1.240 0.215019
## variableASV0040 0.47978 0.69209 0.693 0.488166
## variableASV0042 2.09878 0.54754 3.833 0.000127 ***
## variableASV0044 1.26503 0.59555 2.124 0.033659 *
## variableASV0051 -0.50545 0.50943 -0.992 0.321103
## variableASV0054 0.49609 0.71855 0.690 0.489947
## variableASV0057 0.14685 0.61464 0.239 0.811161
## variableASV0059 0.68076 0.58675 1.160 0.245959
## variableASV0063 0.26277 0.54803 0.479 0.631592
## variableASV0065 -0.30281 0.61917 -0.489 0.624801
## variableASV0066 1.51543 0.72280 2.097 0.036030 *
## variableASV0069 -1.69629 0.62018 -2.735 0.006235 **
## variableASV0072 -0.43166 0.59680 -0.723 0.469506
## variableASV0073 1.70443 0.67058 2.542 0.011031 *
## variableASV0079 -0.49152 0.78820 -0.624 0.532894
## variableASV0080 1.69708 0.69186 2.453 0.014170 *
## variableASV0082 1.15666 0.63788 1.813 0.069785 .
## variableASV0083 -1.25936 0.81448 -1.546 0.122054
## variableASV0084 0.48664 0.62762 0.775 0.438117
## variableASV0085 -0.18806 0.53040 -0.355 0.722916
## variableASV0089 1.01164 0.65379 1.547 0.121782
## variableASV0091 1.06317 0.80740 1.317 0.187911
## variableASV0093 -1.40509 0.60456 -2.324 0.020117 *
## variableASV0094 -0.42514 0.60416 -0.704 0.481628
## variableASV0095 0.50560 0.63661 0.794 0.427074
## variableASV0101 1.64515 0.80063 2.055 0.039897 *
## variableASV0102 0.34004 0.68938 0.493 0.621834
## variableASV0105 1.59824 1.07869 1.482 0.138431
## variableASV0111 -0.95686 0.57235 -1.672 0.094560 .
## variableASV0112 1.24346 0.76621 1.623 0.104617
## variableASV0115 0.53222 0.64288 0.828 0.407742
## variableASV0116 0.03398 0.57100 0.060 0.952550
## variableASV0119 -0.03823 0.62799 -0.061 0.951459
## variableASV0123 0.21105 0.65198 0.324 0.746158
## variableASV0128 1.26038 0.83501 1.509 0.131194
## variableASV0139 1.53228 0.82633 1.854 0.063694 .
## variableASV0145 -1.86479 0.94412 -1.975 0.048250 *
## variableASV0146 0.66931 0.80331 0.833 0.404737
## variableASV0152 1.39689 1.10004 1.270 0.204138
## variableASV0155 -1.17553 0.71711 -1.639 0.101157
## variableASV0167 0.31581 0.69560 0.454 0.649816
## variableASV0175 1.43254 1.10308 1.299 0.194058
## variableASV0180 -1.07281 0.62117 -1.727 0.084152 .
## variableASV0184 1.09999 0.98249 1.120 0.262889
## Microbe_treatmentECO -1.33148 0.52429 -2.540 0.011098 *
## Microbe_treatmentEMO -0.59136 0.14822 -3.990 6.61e-05 ***
## Microbe_treatmentMEM -0.66870 0.14362 -4.656 3.22e-06 ***
## Microbe_treatmentMMO 1.35596 0.24013 5.647 1.63e-08 ***
## Food_typeLarval food -0.29733 0.13983 -2.126 0.033470 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1water.mod.nb2 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom2,dispformula=~variable+Microbe_treatment+Food_type,data=df.all.water.long)
summary(water.mod.nb2)## Family: nbinom2 ( log )
## Formula:
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type)
## Dispersion: ~variable + Microbe_treatment + Food_type
## Data: df.all.water.long
##
## AIC BIC logLik deviance df.resid
## 23799.3 24340.2 -11811.7 23623.3 3362
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 3.444e-06 0.001856
## variable:Microbe_treatment (Intercept) 1.963e+01 4.431121
## variable:Food_type (Intercept) 2.196e+01 4.686453
## variable:Microbe_treatment:Food_type (Intercept) 3.387e-01 0.582020
## Number of obs: 3450, groups:
## variable, 69; variable:Microbe_treatment, 345; variable:Food_type, 138; variable:Microbe_treatment:Food_type, 690
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.70990 0.80867 -8.297 < 2e-16
## Microbe_treatmentECO -13.46838 0.98961 -13.610 < 2e-16
## Microbe_treatmentEMO -0.07208 0.80083 -0.090 0.928285
## Microbe_treatmentMEM -0.25176 0.79989 -0.315 0.752957
## Microbe_treatmentMMO -6.68505 1.22188 -5.471 4.47e-08
## Food_typeLarval food -4.75081 0.91058 -5.217 1.82e-07
## Microbe_treatmentECO:Food_typeLarval food -1.41242 0.78679 -1.795 0.072627
## Microbe_treatmentEMO:Food_typeLarval food 1.24419 0.37240 3.341 0.000835
## Microbe_treatmentMEM:Food_typeLarval food 1.51652 0.37018 4.097 4.19e-05
## Microbe_treatmentMMO:Food_typeLarval food 1.19207 0.60428 1.973 0.048528
##
## (Intercept) ***
## Microbe_treatmentECO ***
## Microbe_treatmentEMO
## Microbe_treatmentMEM
## Microbe_treatmentMMO ***
## Food_typeLarval food ***
## Microbe_treatmentECO:Food_typeLarval food .
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food ***
## Microbe_treatmentMMO:Food_typeLarval food *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.59171 0.35208 -4.521 6.16e-06 ***
## variableASV0002 2.98368 0.40931 7.289 3.11e-13 ***
## variableASV0003 2.29908 0.45337 5.071 3.95e-07 ***
## variableASV0004 2.69638 0.42981 6.273 3.53e-10 ***
## variableASV0005 3.09013 0.43441 7.113 1.13e-12 ***
## variableASV0007 1.94350 0.41292 4.707 2.52e-06 ***
## variableASV0008 1.44625 0.42608 3.394 0.000688 ***
## variableASV0009 3.61752 0.40976 8.828 < 2e-16 ***
## variableASV0010 2.95017 0.43147 6.837 8.06e-12 ***
## variableASV0011 1.53679 0.42889 3.583 0.000339 ***
## variableASV0012 1.67829 0.41171 4.076 4.57e-05 ***
## variableASV0015 0.32703 0.41830 0.782 0.434336
## variableASV0016 4.41497 0.43980 10.039 < 2e-16 ***
## variableASV0017 1.59587 0.49810 3.204 0.001356 **
## variableASV0018 0.20353 0.61338 0.332 0.740028
## variableASV0020 1.26221 0.52218 2.417 0.015640 *
## variableASV0021 -0.19304 0.53055 -0.364 0.715974
## variableASV0022 1.59617 0.42210 3.781 0.000156 ***
## variableASV0023 0.94622 0.51645 1.832 0.066928 .
## variableASV0024 1.55626 0.42474 3.664 0.000248 ***
## variableASV0025 3.07164 0.52663 5.833 5.46e-09 ***
## variableASV0027 2.27425 0.54024 4.210 2.56e-05 ***
## variableASV0033 -1.70585 0.55628 -3.067 0.002165 **
## variableASV0036 -0.89001 0.44272 -2.010 0.044398 *
## variableASV0037 -2.02357 0.49563 -4.083 4.45e-05 ***
## variableASV0039 0.16350 0.54971 0.297 0.766140
## variableASV0040 0.24090 0.52449 0.459 0.646016
## variableASV0042 0.70825 0.41900 1.690 0.090962 .
## variableASV0044 1.97892 0.50990 3.881 0.000104 ***
## variableASV0051 1.34829 0.42171 3.197 0.001388 **
## variableASV0054 -0.79006 0.49229 -1.605 0.108520
## variableASV0057 0.08184 0.48140 0.170 0.865012
## variableASV0059 -0.03230 0.45870 -0.070 0.943862
## variableASV0063 0.51779 0.44137 1.173 0.240743
## variableASV0065 0.62284 0.49265 1.264 0.206137
## variableASV0066 -0.77537 0.49398 -1.570 0.116499
## variableASV0069 3.97236 0.53702 7.397 1.39e-13 ***
## variableASV0072 2.89809 0.52554 5.515 3.50e-08 ***
## variableASV0073 -0.81592 0.47620 -1.713 0.086637 .
## variableASV0079 0.29782 0.55622 0.535 0.592347
## variableASV0080 -0.84527 0.47211 -1.790 0.073391 .
## variableASV0082 -0.52171 0.45860 -1.138 0.255277
## variableASV0083 1.51843 0.59704 2.543 0.010982 *
## variableASV0084 0.40152 0.47282 0.849 0.395763
## variableASV0085 0.94584 0.43706 2.164 0.030455 *
## variableASV0089 0.58796 0.51865 1.134 0.256947
## variableASV0091 -0.82137 0.56218 -1.461 0.144006
## variableASV0093 2.76554 0.52528 5.265 1.40e-07 ***
## variableASV0094 3.42511 0.51985 6.589 4.44e-11 ***
## variableASV0095 0.90323 0.51196 1.764 0.077688 .
## variableASV0101 -1.21368 0.48963 -2.479 0.013184 *
## variableASV0102 0.73940 0.51083 1.447 0.147773
## variableASV0105 -1.04665 0.74319 -1.408 0.159036
## variableASV0111 0.08368 0.46172 0.181 0.856185
## variableASV0112 -0.23926 0.52980 -0.452 0.651548
## variableASV0115 -0.48014 0.49760 -0.965 0.334588
## variableASV0116 -0.11382 0.43731 -0.260 0.794664
## variableASV0119 1.15229 0.52336 2.202 0.027686 *
## variableASV0123 0.64408 0.52161 1.235 0.216909
## variableASV0128 -1.33402 0.52480 -2.542 0.011024 *
## variableASV0139 -0.57802 0.54388 -1.063 0.287886
## variableASV0145 1.77656 0.63592 2.794 0.005211 **
## variableASV0146 -0.72418 0.54123 -1.338 0.180886
## variableASV0152 -1.40699 0.60350 -2.331 0.019733 *
## variableASV0155 3.57915 0.59920 5.973 2.33e-09 ***
## variableASV0167 0.36274 0.52882 0.686 0.492749
## variableASV0175 -1.40188 0.60345 -2.323 0.020173 *
## variableASV0180 1.72267 0.54713 3.149 0.001641 **
## variableASV0184 -1.83812 0.55782 -3.295 0.000984 ***
## Microbe_treatmentECO 3.59384 0.42409 8.474 < 2e-16 ***
## Microbe_treatmentEMO 0.60138 0.12212 4.924 8.46e-07 ***
## Microbe_treatmentMEM 0.49914 0.11656 4.282 1.85e-05 ***
## Microbe_treatmentMMO -1.04862 0.23045 -4.550 5.36e-06 ***
## Food_typeLarval food -0.50118 0.12017 -4.171 3.04e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## water.mod.nb1 0.0 88
## water.mod.nb2 801.1 88
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##dispersal formula checks
water.mod.nb1.disnovar <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~Microbe_treatment+Food_type,data=df.all.water.long)
anova(water.mod.nb1.disnovar,water.mod.nb1) #sig***## Data: df.all.water.long
## Models:
## water.mod.nb1.disnovar: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1.disnovar: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.disnovar: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~Microbe_treatment + Food_type
## water.mod.nb1: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## water.mod.nb1.disnovar 20 23320 23443 -11640 23280
## water.mod.nb1 88 22998 23539 -11411 22822 458.29 68 < 2.2e-16
##
## water.mod.nb1.disnovar
## water.mod.nb1 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##removing 3 way interaction
water.mod.nb1.nofm <- update(water.mod.nb1,.~.-(1|variable:Microbe_treatment:Food_type))
anova(water.mod.nb1,water.mod.nb1.nofm)## Data: df.all.water.long
## Models:
## water.mod.nb1.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.nofm: variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1.nofm: Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb1: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## water.mod.nb1.nofm 87 23065 23600 -11446 22891
## water.mod.nb1 88 22998 23539 -11411 22822 69.225 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##super sig 0.000242
water.mod.nb2.nofm <- update(water.mod.nb2,.~.-(1|variable:Microbe_treatment:Food_type))
anova(water.mod.nb2,water.mod.nb2.nofm)## Data: df.all.water.long
## Models:
## water.mod.nb2.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2.nofm: variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2.nofm: Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Microbe_treatment + Food_type
## water.mod.nb2: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Microbe_treatment + Food_type
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## water.mod.nb2.nofm 87 23825 24360 -11826 23651
## water.mod.nb2 88 23799 24340 -11812 23623 27.894 1 1.281e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##also super sig
###save conditional modes
water.mod.ranef <- as.data.frame(ranef(water.mod.nb1))
#write.csv(water.mod.ranef, "water.mod.ranef.all.csv")
#saveRDS(water.mod.nb1,file="water.mod.nb1.all.rds")3.3.1 Nbinom1 vs. nbinom2 things
Following this wonderful resource, on page 16
Ben Bolker, Mollie Brooks, Beth Gardner, Cleridy Lennert, Mihoko Minami, October 23, 2012, Owls example: a zero-inflated, generalized linear mixed model for count data. https://groups.nceas.ucsb.edu/non-linear-modeling/projects/owls/WRITEUP/owls.pdf
mvtab <- ddply(df.all.water.long,
.(variable:Microbe_treatment:Food_type),
summarise,
callmean=mean(value),
callvar=var(value))
q1 <- qplot(callmean,callvar,data=mvtab)
print(q1+
## linear (quasi-Poisson/NB1) fit
geom_smooth(method="lm",formula=y~x-1,colour="pink")+
## smooth (loess)
geom_smooth(colour="red")+
## semi-quadratic (NB2/LNP)
geom_smooth(method="lm",formula=y~I(x^2)+offset(x)-1,colour="purple")+
## Poisson (v=m)
geom_abline(a=0,b=1,lty=2))## Warning in geom_abline(a = 0, b = 1, lty = 2): Ignoring unknown parameters: `a`
## and `b`## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
3.4 Larvae - all treatments
## 'data.frame': 2700 obs. of 8 variables:
## $ sum : num 63934 18777 13503 111 2064 ...
## $ Mesocosm_id : chr "B_EMO.3" "B_EMO.4" "B_EMO.5" "B_ECO.1" ...
## $ Microbe_treatment: chr "EMO" "EMO" "EMO" "ECO" ...
## $ Food_type : chr "Bamboo" "Bamboo" "Bamboo" "Bamboo" ...
## $ Biomass_day40 : num 0.621 0.308 4.085 0 0 ...
## $ variable : Factor w/ 54 levels "ASV0001","ASV0002",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 0 0 15 111 1867 ...
## $ rowid : int 1 2 3 4 5 6 7 8 9 10 ...df.all.larv.long$rowid <- as.factor(df.all.larv.long$rowid)
df.all.larv.long$Microbe_treatment <- as.factor(df.all.larv.long$Microbe_treatment)
df.all.larv.long$Food_type <- as.factor(df.all.larv.long$Food_type)
#df.all.larv.long$food.microbes <- as.factor(df.all.larv.long$food.microbes)
df.all.larv.long$Mesocosm_id <- as.factor(df.all.larv.long$Mesocosm_id)
str(df.all.larv.long)## 'data.frame': 2700 obs. of 8 variables:
## $ sum : num 63934 18777 13503 111 2064 ...
## $ Mesocosm_id : Factor w/ 50 levels "B_ALL.1","B_ALL.2",..: 13 14 15 6 7 8 9 10 16 17 ...
## $ Microbe_treatment: Factor w/ 5 levels "ALL","ECO","EMO",..: 3 3 3 2 2 2 2 2 4 4 ...
## $ Food_type : Factor w/ 2 levels "Bamboo","Larval food": 1 1 1 1 1 1 1 1 1 1 ...
## $ Biomass_day40 : num 0.621 0.308 4.085 0 0 ...
## $ variable : Factor w/ 54 levels "ASV0001","ASV0002",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 0 0 15 111 1867 ...
## $ rowid : Factor w/ 2700 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...##binomial
larv.mod.bin <- glmmTMB(cbind(value, sum-value)~Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type)+(1|rowid),family=binomial,data=df.all.larv.long)
summary(larv.mod.bin)## Family: binomial ( logit )
## Formula:
## cbind(value, sum - value) ~ Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type) + (1 | rowid)
## Data: df.all.larv.long
##
## AIC BIC logLik deviance df.resid
## 15583.1 15671.6 -7776.6 15553.1 2685
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 8.950e-07 0.0009461
## variable:Microbe_treatment (Intercept) 2.604e+01 5.1032341
## variable:Food_type (Intercept) 3.717e+01 6.0969133
## variable:Microbe_treatment:Food_type (Intercept) 6.807e+00 2.6090817
## rowid (Intercept) 1.322e+01 3.6365737
## Number of obs: 2700, groups:
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108; variable:Microbe_treatment:Food_type, 540; rowid, 2700
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.4884 1.2226 -9.397 < 2e-16
## Microbe_treatmentECO -8.9465 1.3737 -6.513 7.37e-11
## Microbe_treatmentEMO -2.8262 1.2533 -2.255 0.024137
## Microbe_treatmentMEM 0.1802 1.2341 0.146 0.883936
## Microbe_treatmentMMO -8.2294 1.2300 -6.691 2.22e-11
## Food_typeLarval food -5.8031 1.4682 -3.953 7.73e-05
## Microbe_treatmentECO:Food_typeLarval food -2.3232 1.8710 -1.242 0.214342
## Microbe_treatmentEMO:Food_typeLarval food 3.8919 1.1002 3.538 0.000404
## Microbe_treatmentMEM:Food_typeLarval food 1.6363 1.0968 1.492 0.135748
## Microbe_treatmentMMO:Food_typeLarval food 2.8462 1.2093 2.354 0.018593
##
## (Intercept) ***
## Microbe_treatmentECO ***
## Microbe_treatmentEMO *
## Microbe_treatmentMEM
## Microbe_treatmentMMO ***
## Food_typeLarval food ***
## Microbe_treatmentECO:Food_typeLarval food
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food
## Microbe_treatmentMMO:Food_typeLarval food *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##full model
larv.mod.nb1 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),disp=~variable+Food_type+Microbe_treatment,family=nbinom1,data=df.all.larv.long)
summary(larv.mod.nb1)## Family: nbinom1 ( log )
## Formula:
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type)
## Dispersion: ~variable + Food_type + Microbe_treatment
## Data: df.all.larv.long
##
## AIC BIC logLik deviance df.resid
## 14676.1 15106.9 -7265.0 14530.1 2627
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 1.860e-07 0.0004313
## variable:Microbe_treatment (Intercept) 3.833e+00 1.9577005
## variable:Food_type (Intercept) 6.559e+00 2.5610822
## variable:Microbe_treatment:Food_type (Intercept) 2.367e-01 0.4864757
## Number of obs: 2700, groups:
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108; variable:Microbe_treatment:Food_type, 540
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.4933 0.4805 -13.513 < 2e-16
## Microbe_treatmentECO -2.8758 0.9653 -2.979 0.002889
## Microbe_treatmentEMO -0.5272 0.4613 -1.143 0.253106
## Microbe_treatmentMEM -0.1958 0.4377 -0.447 0.654578
## Microbe_treatmentMMO -2.1999 0.5132 -4.287 1.81e-05
## Food_typeLarval food -1.1961 0.5968 -2.004 0.045035
## Microbe_treatmentECO:Food_typeLarval food -2.3653 0.6487 -3.646 0.000266
## Microbe_treatmentEMO:Food_typeLarval food 0.9873 0.3431 2.877 0.004010
## Microbe_treatmentMEM:Food_typeLarval food 0.7123 0.3163 2.252 0.024326
## Microbe_treatmentMMO:Food_typeLarval food 1.1627 0.3863 3.010 0.002611
##
## (Intercept) ***
## Microbe_treatmentECO **
## Microbe_treatmentEMO
## Microbe_treatmentMEM
## Microbe_treatmentMMO ***
## Food_typeLarval food *
## Microbe_treatmentECO:Food_typeLarval food ***
## Microbe_treatmentEMO:Food_typeLarval food **
## Microbe_treatmentMEM:Food_typeLarval food *
## Microbe_treatmentMMO:Food_typeLarval food **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.608016 0.480928 13.740 < 2e-16 ***
## variableASV0002 0.125802 0.830243 0.152 0.879562
## variableASV0003 1.424233 0.541764 2.629 0.008567 **
## variableASV0004 1.377207 0.989981 1.391 0.164181
## variableASV0005 0.568316 0.729597 0.779 0.436012
## variableASV0006 1.481056 0.556652 2.661 0.007799 **
## variableASV0007 0.205723 0.574432 0.358 0.720244
## variableASV0009 0.805902 0.583140 1.382 0.166971
## variableASV0010 0.120688 0.643408 0.188 0.851209
## variableASV0011 -0.744663 0.627092 -1.187 0.235036
## variableASV0012 0.397982 0.814716 0.488 0.625202
## variableASV0013 1.423164 0.668655 2.128 0.033304 *
## variableASV0016 0.710856 0.765013 0.929 0.352781
## variableASV0019 1.814936 0.676409 2.683 0.007292 **
## variableASV0024 -0.120270 0.592496 -0.203 0.839144
## variableASV0027 0.616208 0.640521 0.962 0.336029
## variableASV0030 -0.101724 1.076445 -0.095 0.924712
## variableASV0032 0.452840 0.710405 0.637 0.523839
## variableASV0033 1.277173 0.603836 2.115 0.034422 *
## variableASV0043 3.361423 0.847375 3.967 7.28e-05 ***
## variableASV0053 0.366468 0.629382 0.582 0.560387
## variableASV0054 0.409139 0.820909 0.498 0.618204
## variableASV0055 1.455563 0.696111 2.091 0.036529 *
## variableASV0058 1.669642 0.711720 2.346 0.018980 *
## variableASV0059 -0.524573 0.607181 -0.864 0.387616
## variableASV0060 2.074881 0.759824 2.731 0.006319 **
## variableASV0061 1.279328 0.660277 1.938 0.052677 .
## variableASV0066 -1.241301 0.686333 -1.809 0.070513 .
## variableASV0070 -0.060359 0.812150 -0.074 0.940756
## variableASV0079 -0.709432 0.587382 -1.208 0.227129
## variableASV0083 -0.604105 0.588541 -1.026 0.304681
## variableASV0085 0.009405 0.958822 0.010 0.992174
## variableASV0087 0.545854 0.655877 0.832 0.405268
## variableASV0099 0.455684 0.725616 0.628 0.530006
## variableASV0100 -0.289356 0.659121 -0.439 0.660660
## variableASV0102 -1.381170 0.862920 -1.601 0.109471
## variableASV0109 0.596347 1.010469 0.590 0.555078
## variableASV0125 1.346848 0.776236 1.735 0.082723 .
## variableASV0130 0.076914 0.678606 0.113 0.909760
## variableASV0134 0.931249 0.904047 1.030 0.302968
## variableASV0138 -0.864209 0.741853 -1.165 0.244046
## variableASV0149 -2.035605 0.624026 -3.262 0.001106 **
## variableASV0150 -1.340375 0.689884 -1.943 0.052028 .
## variableASV0152 0.891285 0.907213 0.982 0.325882
## variableASV0164 -1.121652 0.661599 -1.695 0.090006 .
## variableASV0210 0.365497 1.072546 0.341 0.733273
## variableASV0212 0.620222 0.959904 0.646 0.518196
## variableASV0218 -1.112410 0.845764 -1.315 0.188418
## variableASV0220 1.082720 1.052458 1.029 0.303595
## variableASV0225 0.297385 1.068749 0.278 0.780816
## variableASV0262 -0.372923 1.053183 -0.354 0.723271
## variableASV0276 0.684188 1.024039 0.668 0.504053
## variableASV0289 -1.178319 0.801413 -1.470 0.141480
## variableASV0311 -0.957374 0.776274 -1.233 0.217466
## Food_typeLarval food 1.045981 0.222090 4.710 2.48e-06 ***
## Microbe_treatmentECO -1.008090 0.807555 -1.248 0.211913
## Microbe_treatmentEMO -0.100263 0.223110 -0.449 0.653152
## Microbe_treatmentMEM -0.661644 0.194947 -3.394 0.000689 ***
## Microbe_treatmentMMO 0.208428 0.263089 0.792 0.428224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1larv.mod.nb2 <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),disp=~variable+Food_type+Microbe_treatment,family=nbinom2,data=df.all.larv.long)
summary(larv.mod.nb2)## Family: nbinom2 ( log )
## Formula:
## value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | variable:Microbe_treatment:Food_type)
## Dispersion: ~variable + Food_type + Microbe_treatment
## Data: df.all.larv.long
##
## AIC BIC logLik deviance df.resid
## 15247.7 15678.5 -7550.9 15101.7 2627
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 1.099e-06 0.001048
## variable:Microbe_treatment (Intercept) 1.389e+01 3.727191
## variable:Food_type (Intercept) 2.584e+01 5.082913
## variable:Microbe_treatment:Food_type (Intercept) 1.727e+00 1.313982
## Number of obs: 2700, groups:
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108; variable:Microbe_treatment:Food_type, 540
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.2123 0.9771 -5.334 9.58e-08
## Microbe_treatmentECO -10.3785 1.1809 -8.789 < 2e-16
## Microbe_treatmentEMO -3.4299 0.9398 -3.649 0.000263
## Microbe_treatmentMEM -1.0148 0.9043 -1.122 0.261784
## Microbe_treatmentMMO -5.4497 1.0973 -4.967 6.81e-07
## Food_typeLarval food -5.6533 1.1879 -4.759 1.95e-06
## Microbe_treatmentECO:Food_typeLarval food -1.2564 1.3606 -0.923 0.355788
## Microbe_treatmentEMO:Food_typeLarval food 4.3010 0.8108 5.305 1.13e-07
## Microbe_treatmentMEM:Food_typeLarval food 1.6533 0.7612 2.172 0.029853
## Microbe_treatmentMMO:Food_typeLarval food 1.8631 0.9013 2.067 0.038718
##
## (Intercept) ***
## Microbe_treatmentECO ***
## Microbe_treatmentEMO ***
## Microbe_treatmentMEM
## Microbe_treatmentMMO ***
## Food_typeLarval food ***
## Microbe_treatmentECO:Food_typeLarval food
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food *
## Microbe_treatmentMMO:Food_typeLarval food *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.243182 0.350061 -9.265 < 2e-16 ***
## variableASV0002 -0.591657 0.452883 -1.306 0.191408
## variableASV0003 1.985713 0.401882 4.941 7.77e-07 ***
## variableASV0004 -1.074704 0.554580 -1.938 0.052639 .
## variableASV0005 0.066652 0.508792 0.131 0.895775
## variableASV0006 2.174526 0.438935 4.954 7.27e-07 ***
## variableASV0007 1.385602 0.440130 3.148 0.001643 **
## variableASV0009 0.969397 0.432365 2.242 0.024956 *
## variableASV0010 0.007023 0.450935 0.016 0.987574
## variableASV0011 1.343682 0.501763 2.678 0.007408 **
## variableASV0012 -0.493686 0.508836 -0.970 0.331934
## variableASV0013 4.735712 0.524659 9.026 < 2e-16 ***
## variableASV0016 -0.193298 0.515311 -0.375 0.707579
## variableASV0019 0.494949 0.496496 0.997 0.318820
## variableASV0024 1.352986 0.436170 3.102 0.001922 **
## variableASV0027 2.222017 0.485900 4.573 4.81e-06 ***
## variableASV0030 -2.229292 0.576144 -3.869 0.000109 ***
## variableASV0032 4.293441 0.517095 8.303 < 2e-16 ***
## variableASV0033 0.846227 0.411874 2.055 0.039920 *
## variableASV0043 0.310331 0.567328 0.547 0.584375
## variableASV0053 1.535773 0.463201 3.316 0.000915 ***
## variableASV0054 -0.727097 0.478087 -1.521 0.128298
## variableASV0055 2.429137 0.489142 4.966 6.83e-07 ***
## variableASV0058 2.095440 0.500062 4.190 2.79e-05 ***
## variableASV0059 1.057266 0.450351 2.348 0.018892 *
## variableASV0060 0.889673 0.489455 1.818 0.069113 .
## variableASV0061 0.565403 0.438701 1.289 0.197463
## variableASV0066 0.250592 0.500883 0.500 0.616864
## variableASV0070 2.241231 0.536312 4.179 2.93e-05 ***
## variableASV0079 2.149780 0.425556 5.052 4.38e-07 ***
## variableASV0083 1.649553 0.423276 3.897 9.73e-05 ***
## variableASV0085 -0.793906 0.578679 -1.372 0.170086
## variableASV0087 3.098546 0.494511 6.266 3.71e-10 ***
## variableASV0099 -0.147528 0.446542 -0.330 0.741114
## variableASV0100 1.045376 0.506538 2.064 0.039040 *
## variableASV0102 -0.489850 0.572942 -0.855 0.392565
## variableASV0109 0.651864 0.620928 1.050 0.293799
## variableASV0125 1.137772 0.495724 2.295 0.021723 *
## variableASV0130 0.213252 0.445147 0.479 0.631897
## variableASV0134 -0.683177 0.520415 -1.313 0.189266
## variableASV0138 1.987035 0.525266 3.783 0.000155 ***
## variableASV0149 3.648389 0.495782 7.359 1.85e-13 ***
## variableASV0150 1.336130 0.512202 2.609 0.009091 **
## variableASV0152 0.804509 0.569131 1.414 0.157486
## variableASV0164 0.887599 0.477293 1.860 0.062934 .
## variableASV0210 -1.296908 0.576441 -2.250 0.024458 *
## variableASV0212 0.758784 0.591754 1.282 0.199751
## variableASV0218 0.115050 0.581424 0.198 0.843141
## variableASV0220 0.392096 0.625029 0.627 0.530446
## variableASV0225 0.452762 0.647796 0.699 0.484598
## variableASV0262 0.274559 0.787801 0.349 0.727455
## variableASV0276 0.343902 0.625959 0.549 0.582730
## variableASV0289 1.540508 0.559074 2.755 0.005861 **
## variableASV0311 1.767208 0.540069 3.272 0.001067 **
## Food_typeLarval food 0.967602 0.130601 7.409 1.27e-13 ***
## Microbe_treatmentECO 2.851151 0.397070 7.180 6.95e-13 ***
## Microbe_treatmentEMO 0.355197 0.147717 2.405 0.016191 *
## Microbe_treatmentMEM 0.523868 0.136667 3.833 0.000127 ***
## Microbe_treatmentMMO -0.143964 0.205901 -0.699 0.484434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## larv.mod.nb1 0.0 73
## larv.mod.nb2 571.6 73
## larv.mod.bin 907.0 15
##dispersal formula checks
larv.mod.nb1.disnovar <- glmmTMB(value~offset(log(sum))+Microbe_treatment*Food_type+(1|variable)+(1|variable:Microbe_treatment)+(1|variable:Food_type)+(1|variable:Microbe_treatment:Food_type),family=nbinom1,dispformula=~Microbe_treatment+Food_type,data=df.all.larv.long)
anova(larv.mod.nb1.disnovar,larv.mod.nb1) #sig***## Data: df.all.larv.long
## Models:
## larv.mod.nb1.disnovar: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1.disnovar: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.disnovar: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~Microbe_treatment + Food_type
## larv.mod.nb1: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## larv.mod.nb1.disnovar 20 14775 14893 -7367.4 14735
## larv.mod.nb1 73 14676 15107 -7265.0 14530 204.7 53 < 2.2e-16
##
## larv.mod.nb1.disnovar
## larv.mod.nb1 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##removing 3 way interaction
larv.mod.nb1.nofm <- update(larv.mod.nb1,.~.-(1|variable:Microbe_treatment:Food_type))
summary(larv.mod.nb1.nofm)## Family: nbinom1 ( log )
## Formula:
## value ~ Microbe_treatment + Food_type + (1 | variable) + (1 |
## variable:Microbe_treatment) + (1 | variable:Food_type) +
## Microbe_treatment:Food_type + offset(log(sum))
## Dispersion: ~variable + Food_type + Microbe_treatment
## Data: df.all.larv.long
##
## AIC BIC logLik deviance df.resid
## 14678.5 15103.4 -7267.3 14534.5 2628
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 1.312e-07 0.0003623
## variable:Microbe_treatment (Intercept) 4.029e+00 2.0072916
## variable:Food_type (Intercept) 6.676e+00 2.5837467
## Number of obs: 2700, groups:
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.4678 0.4802 -13.469 < 2e-16
## Microbe_treatmentECO -2.8162 0.9564 -2.945 0.003233
## Microbe_treatmentEMO -0.5032 0.4541 -1.108 0.267817
## Microbe_treatmentMEM -0.1419 0.4287 -0.331 0.740638
## Microbe_treatmentMMO -2.2053 0.4972 -4.435 9.19e-06
## Food_typeLarval food -1.2519 0.5825 -2.149 0.031613
## Microbe_treatmentECO:Food_typeLarval food -2.6340 0.5224 -5.042 4.61e-07
## Microbe_treatmentEMO:Food_typeLarval food 0.9767 0.2841 3.438 0.000586
## Microbe_treatmentMEM:Food_typeLarval food 0.7184 0.2366 3.036 0.002395
## Microbe_treatmentMMO:Food_typeLarval food 1.1724 0.2773 4.228 2.36e-05
##
## (Intercept) ***
## Microbe_treatmentECO **
## Microbe_treatmentEMO
## Microbe_treatmentMEM
## Microbe_treatmentMMO ***
## Food_typeLarval food *
## Microbe_treatmentECO:Food_typeLarval food ***
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food **
## Microbe_treatmentMMO:Food_typeLarval food ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.85166 0.46809 14.638 < 2e-16 ***
## variableASV0002 0.04264 0.85391 0.050 0.960178
## variableASV0003 1.17327 0.53168 2.207 0.027334 *
## variableASV0004 1.15220 0.98560 1.169 0.242390
## variableASV0005 0.38587 0.72643 0.531 0.595292
## variableASV0006 1.35183 0.56042 2.412 0.015857 *
## variableASV0007 0.02670 0.57158 0.047 0.962747
## variableASV0009 0.61519 0.57860 1.063 0.287675
## variableASV0010 -0.07145 0.63805 -0.112 0.910843
## variableASV0011 -0.92455 0.62355 -1.483 0.138146
## variableASV0012 0.20960 0.81152 0.258 0.796190
## variableASV0013 1.14764 0.65803 1.744 0.081149 .
## variableASV0016 0.52626 0.76301 0.690 0.490370
## variableASV0019 1.58854 0.66730 2.381 0.017287 *
## variableASV0024 -0.33074 0.58623 -0.564 0.572629
## variableASV0027 0.36266 0.63071 0.575 0.565290
## variableASV0030 -0.14127 1.08752 -0.130 0.896644
## variableASV0032 0.21191 0.70540 0.300 0.763861
## variableASV0033 1.08423 0.60075 1.805 0.071105 .
## variableASV0043 3.09257 0.84087 3.678 0.000235 ***
## variableASV0053 0.10922 0.61851 0.177 0.859838
## variableASV0054 0.26393 0.83115 0.318 0.750833
## variableASV0055 1.21154 0.68891 1.759 0.078641 .
## variableASV0058 1.41797 0.70361 2.015 0.043875 *
## variableASV0059 -0.71776 0.60315 -1.190 0.234036
## variableASV0060 1.82663 0.75287 2.426 0.015256 *
## variableASV0061 1.03000 0.64992 1.585 0.113010
## variableASV0066 -1.42025 0.68413 -2.076 0.037896 *
## variableASV0070 -0.27711 0.81396 -0.340 0.733523
## variableASV0079 -0.95055 0.57573 -1.651 0.098736 .
## variableASV0083 -0.77493 0.58019 -1.336 0.181665
## variableASV0085 -0.18678 0.95542 -0.195 0.845007
## variableASV0087 0.27474 0.64468 0.426 0.669991
## variableASV0099 0.24630 0.72144 0.341 0.732797
## variableASV0100 -0.47850 0.65498 -0.731 0.465046
## variableASV0102 -1.57001 0.85924 -1.827 0.067667 .
## variableASV0109 0.32935 1.00606 0.327 0.743391
## variableASV0125 1.08067 0.76783 1.407 0.159302
## variableASV0130 -0.01441 0.68679 -0.021 0.983255
## variableASV0134 0.73435 0.90294 0.813 0.416053
## variableASV0138 -1.13267 0.73269 -1.546 0.122123
## variableASV0149 -2.29236 0.61342 -3.737 0.000186 ***
## variableASV0150 -1.53318 0.68541 -2.237 0.025293 *
## variableASV0152 0.61808 0.89867 0.688 0.491598
## variableASV0164 -1.33252 0.65650 -2.030 0.042382 *
## variableASV0210 0.17882 1.06693 0.168 0.866893
## variableASV0212 0.35378 0.95374 0.371 0.710678
## variableASV0218 -1.30772 0.84022 -1.556 0.119611
## variableASV0220 0.81224 1.04553 0.777 0.437237
## variableASV0225 0.02016 1.06207 0.019 0.984859
## variableASV0262 -0.56844 1.04866 -0.542 0.587775
## variableASV0276 0.43079 1.02047 0.422 0.672915
## variableASV0289 -1.44067 0.79490 -1.812 0.069925 .
## variableASV0311 -1.22539 0.76813 -1.595 0.110647
## Food_typeLarval food 0.96112 0.21909 4.387 1.15e-05 ***
## Microbe_treatmentECO -1.07382 0.82793 -1.297 0.194634
## Microbe_treatmentEMO -0.07345 0.22125 -0.332 0.739912
## Microbe_treatmentMEM -0.61978 0.19197 -3.229 0.001244 **
## Microbe_treatmentMMO 0.21159 0.26152 0.809 0.418463
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.all.larv.long
## Models:
## larv.mod.nb1.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.nofm: variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1.nofm: Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb1: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## larv.mod.nb1.nofm 72 14678 15103 -7267.3 14534
## larv.mod.nb1 73 14676 15107 -7265.0 14530 4.4686 1 0.03452 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ##removing just food
# larv.mod.nb1.nofood <- update(larv.mod.nb1,.~.-(1|variable:Food_type))
# anova(larv.mod.nb1,larv.mod.nb1.nofood) #super sig***
#
# ##removing just microbes
# larv.mod.nb1.nomicr <- update(larv.mod.nb1,.~.-(1|variable:Microbe_treatment))
# anova(larv.mod.nb1,larv.mod.nb1.nomicr) #super sig***
##checking out nbinom2
larv.mod.nb2.nofm <- update(larv.mod.nb2,.~.-(1|variable:Microbe_treatment:Food_type))
anova(larv.mod.nb2,larv.mod.nb2.nofm)## Data: df.all.larv.long
## Models:
## larv.mod.nb2.nofm: value ~ Microbe_treatment + Food_type + (1 | variable) + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2.nofm: variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2.nofm: Microbe_treatment:Food_type + offset(log(sum)), zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2: value ~ offset(log(sum)) + Microbe_treatment * Food_type + (1 | , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~variable + Food_type + Microbe_treatment
## larv.mod.nb2: (1 | variable:Microbe_treatment:Food_type), zi=~0, disp=~variable + Food_type + Microbe_treatment
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## larv.mod.nb2.nofm 72 15256 15681 -7556.1 15112
## larv.mod.nb2 73 15248 15678 -7550.9 15102 10.524 1 0.001179 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##sig**
##checking out binomial
larv.mod.bin.nofm <- update(larv.mod.bin,.~.-(1|variable:Microbe_treatment:Food_type))
summary(larv.mod.bin.nofm)## Family: binomial ( logit )
## Formula:
## cbind(value, sum - value) ~ Microbe_treatment + Food_type + (1 |
## variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) +
## (1 | rowid) + Microbe_treatment:Food_type
## Data: df.all.larv.long
##
## AIC BIC logLik deviance df.resid
## 15601.7 15684.3 -7786.9 15573.7 2686
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## variable (Intercept) 1.307e-06 0.001143
## variable:Microbe_treatment (Intercept) 3.089e+01 5.557434
## variable:Food_type (Intercept) 4.042e+01 6.357299
## rowid (Intercept) 1.406e+01 3.749305
## Number of obs: 2700, groups:
## variable, 54; variable:Microbe_treatment, 270; variable:Food_type, 108; rowid, 2700
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.3323 1.2193 -9.294 < 2e-16
## Microbe_treatmentECO -9.0397 1.3069 -6.917 4.62e-12
## Microbe_treatmentEMO -2.6017 1.2026 -2.163 0.03051
## Microbe_treatmentMEM 0.4579 1.1676 0.392 0.69494
## Microbe_treatmentMMO -8.3851 1.1753 -7.135 9.70e-13
## Food_typeLarval food -5.8334 1.4727 -3.961 7.46e-05
## Microbe_treatmentECO:Food_typeLarval food -2.8321 1.4541 -1.948 0.05146
## Microbe_treatmentEMO:Food_typeLarval food 3.5878 0.7640 4.696 2.65e-06
## Microbe_treatmentMEM:Food_typeLarval food 1.8169 0.7501 2.422 0.01543
## Microbe_treatmentMMO:Food_typeLarval food 2.9328 0.8917 3.289 0.00101
##
## (Intercept) ***
## Microbe_treatmentECO ***
## Microbe_treatmentEMO *
## Microbe_treatmentMEM
## Microbe_treatmentMMO ***
## Food_typeLarval food ***
## Microbe_treatmentECO:Food_typeLarval food .
## Microbe_treatmentEMO:Food_typeLarval food ***
## Microbe_treatmentMEM:Food_typeLarval food *
## Microbe_treatmentMMO:Food_typeLarval food **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.all.larv.long
## Models:
## larv.mod.bin.nofm: cbind(value, sum - value) ~ Microbe_treatment + Food_type + (1 | , zi=~0, disp=~1
## larv.mod.bin.nofm: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~1
## larv.mod.bin.nofm: (1 | rowid) + Microbe_treatment:Food_type, zi=~0, disp=~1
## larv.mod.bin: cbind(value, sum - value) ~ Microbe_treatment * Food_type + (1 | , zi=~0, disp=~1
## larv.mod.bin: variable) + (1 | variable:Microbe_treatment) + (1 | variable:Food_type) + , zi=~0, disp=~1
## larv.mod.bin: (1 | variable:Microbe_treatment:Food_type) + (1 | rowid), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## larv.mod.bin.nofm 14 15602 15684 -7786.9 15574
## larv.mod.bin 15 15583 15672 -7776.6 15553 20.617 1 5.61e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1###save conditional modes
larv.mod.ranef <- as.data.frame(ranef(larv.mod.nb1))
write.csv(larv.mod.ranef, "larv.mod.ranef.all.csv")
saveRDS(larv.mod.nb1,file="larv.mod.nb1.all.rds")
# library("ggstats")
# ggcoef_multicomponents(larv.mod.nb1)3.4.1 Nbinom1 vs. nbinom2 things
##note: some errors/warnings if the treatments aren't factors for some reason
mvtab <- ddply(df.all.larv.long,
.(variable:Microbe_treatment:Food_type),
summarise,
callmean=mean(value),
callvar=var(value))
q1 <- qplot(callmean,callvar,data=mvtab)
print(q1+
## linear (quasi-Poisson/NB1) fit
geom_smooth(method="lm",formula=y~x-1,colour="pink")+
## smooth (loess)
geom_smooth(colour="red")+
## semi-quadratic (NB2/LNP)
geom_smooth(method="lm",formula=y~I(x^2)+offset(x)-1,colour="purple")+
## Poisson (v=m)
geom_abline(a=0,b=1,lty=2))#+## Warning in geom_abline(a = 0, b = 1, lty = 2): Ignoring unknown parameters: `a`
## and `b`## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'